Interior point method for reformulated optimal power flow model

ABSTRACT

A method for approximating an optimal power flow of a smart electric power grid includes providing a cost function that models a smart electric power grid having buses connected by branches, deriving a set of linear equations that minimize the cost function subject to constraints from an expression of an extremum of the cost function with respect to all arguments, reducing a dimension of the linear equations by solving for a subset of the linear equations, re-organizing the reduced dimension linear equations into primal and dual parts, and decomposing the re-organized reduced dimensional linear equations into two systems of block matrix equations which can be solved by a series of back substitutions. A solution of the two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through the smart electric power grid.

TECHNICAL FIELD

This disclosure is directed to methods of analyzing and modeling electric power grids, in particular smart grids.

DISCUSSION OF THE RELATED ART

An Optimal Power Flow (OPF) model of an electric power network is a model whose solution provides the conditions which give the lowest cost per kilowatthour delivered. Electricity flowing through an AC network is subject to Kirchhoff's Laws, and transmission lines are subject to thermal limits as well as voltage and electrical stability constraints. An OPF model should be capable to quickly and efficiently analyze an electric power network and yield a feasible power flow solution under a wide variety of adverse conditions and determine the most efficient set of actions to take.

OPF models present challenges even for cutting-edge mathematical programming techniques due to a sheer size of contingency cases that fully couple all network components for all adverse scenarios. Decoupling presents issues as it is well known that contingency analysis of power networks lacks even the most basic intuitive properties such as monotonicity, in that if tripping a branch results in an infeasible network, then tripping this branch and an extra branch should also be infeasible, but often is not.

A basecase OPF model seeks to produce an optimal power flow solution balancing generation/consumption for all buses with power distribution across all connected branches such that the overall system state, as measured by underlying physics parameters including voltages and phase angles, is feasible and safe.

A basecase OPF model can control AC and DC energy dispatches, regulate AC bus voltages and phase angles, select AC transformer/phasor tap and shunt switches, and control the AC/DC converter.

In a basecase OPF model, an electric power network is made up of the following components and their controls/states:

-   -   AC Bus: represented by a node in the network graph. An AC bus         has the following controls & states:         -   Load: Utilizes active (MW) and/or reactive (MVar) power;         -   Generator: Provides active (MW) and reactive (MVar) power;         -   Voltage and Phase angle: restricted to certain range around             nominal;         -   Shunt capacitor/inductor: are optional, can be fixed or             switching, can be positive (capacitor) or negative             (inductor) and moderates active and reactive net power             produced at the bus;         -   The main AC bus parameters are: bounds for Load and             Generation, and bounds for Voltage and Phase Angle.     -   AC Branch/Line: represented by an arc in the network graph and         transmits both active and reactive power. An AC branch has the         following controls and states:         -   Power flow: active and reactive power injected into the             branch by both of its connected buses;         -   Adjusted impedance: is optional and represents branch             impedance adjusted by an attached transformer;         -   The main AC branch parameters are: Admittance matrix,             calculated based on the branch Impedance and the branch             fixed shunt admittances, and Flow Capacity Ratings, which             are the minimum and maximum total flow as constrained by             thermal ratings.     -   AC Transformer: is an extra add-on to the AC branch that may         either transform its connected buses' voltages, (referred to         herein as Transformer mode) or shift the difference between its         connected buses' phase angles (referred to herein below as         Phasor mode). The AC Transformer has the following controls and         states:         -   Tap choice;         -   The main AC Transformer parameters are Tap Turns Ratio for             Transformer mode, and Tap Phase Shift for Phasor mode.     -   DC Bus: represented by a node in the network graph. A DC bus has         the following controls & states:         -   DC Voltage and Current Injection;         -   The main DC bus parameters are: bounds for Voltage and             Current.     -   DC Branch/Line: represented by an arc in the network graph and         transmits a constant current. A DC branch has the following         controls and states:         -   DC Current;         -   The main DC bus parameter is: Resistance.     -   Converter: represented by an arc in the network graph that         connects AC and DC buses. Converters can be of two types:         Rectifier, for AC to DC, and Inverter for DC to AC. A converter         has the following controls and states:         -   Firing and Overlap angles, AC side voltage, since the             converter comes with a Transformer;         -   Power Factor, Active Power injection;         -   The main Converter parameters are: bounds for all converter             controls and states and Commutating Impedance.

Other model components include a Voltage Source Converter (VSC), which is a special type of the converter that does not involve firing/overlap angles and directly controls voltage, power factor and active power injection, and FACTS devices, which are not modeled as distinct components but rather define tight bounds on certain controls/states of network components that are regulated by FACTS. Buses may be grouped in areas for defining area power exchange and inter-area power transfer constraints. Area involves several Interfaces that are corresponding groupings of branches that connect one area to another.

FIG. 1 is labelled to show each of the elements that may exist in an electric power system network. Their symbolic representation may also be seen.

An OPF mathematical formulation according to an embodiment of the invention is focused on making a base-case model formulation suitable for optimization solvers, in particular, utilizing as few as possible non-linear equations and as few as possible integer variables, and is extendable to contingency cases.

A power network is modeled as a standard graph N={V(N), E(N)} where V(N) is the set of buses and E(N) is the set of branches. Within V(N), one can distinguish the following subsets: the AC bus set V(AC) and the DC bus set V(DC). Within E(N) one can distinguish the set of AC branches E(AC), the set of DC branches E(DC), and the set of converters V(C). Also, one cane distinguish inter-area interface branch groupings E(I).

Every control and state variable mentioned in the previous section has a box constraint that are represented by lower and upper bound constraints, not all of which are explicitly mentioned below. An example of box constraint is the active power generation control variable P_(ig), which is subject to generation bound box constraints: P _(ig min) ≦P _(ig) ≦P _(ig max). Note that all choice variables in an OPF model according to an embodiment of the invention, which reflect the choice of each transformer/phasor tap or each shunt switch, denoted w_(k), not only are subject to be within the [0; 1] box but also should satisfy the Exclusive Choice constraint:

${\sum\limits_{k = 1}^{K}\; w_{k}} = 1$ The set of AC bus constraints comprises the fundamental set of network constraints, called the Power Flow Conservation set (PFC):

${P_{ig} - P_{il} + P_{i,{shunt}}} = {{\sum\limits_{e \in {\delta^{+}{(i)}}}\; P_{e}} + {\sum\limits_{c \in {\delta^{+}{(i)}}}\; P_{c}}}$ ${Q_{ig} - Q_{il} + Q_{i,{shunt}}} = {{\sum\limits_{e \in {\delta^{+}{(i)}}}\; Q_{e}} + {\sum\limits_{c \in {\delta^{+}{(i)}}}\;{Q_{c}.}}}$ The l.h.s. reflects the net power injection generation and the r.h.s. reflects the net power flow coming from both V(AC) and V(C). The shunt term in the PFC constraint is given by:

$P_{i,{shunt}} = {\sum\limits_{k = 1}^{K_{i}}\;{w_{ik} \cdot P_{ik}}}$ $Q_{i,{shunt}} = {\sum\limits_{k = 1}^{K_{i}}\;{w_{ik} \cdot Q_{ik}}}$ The set of AC branch constraints includes:

${{Voltage}\mspace{14mu}{Law}\mspace{11mu}({VL})}:\begin{matrix} {P_{ij} = {{V_{ii}G_{ii}} + {V_{ij}\left( {{G_{ij}*{\cos\left( \theta_{ij} \right)}} + {B_{ij}*{\sin\left( \theta_{ij} \right)}}} \right)}}} \\ {Q_{ij} = {{{- V_{ii}}B_{ii}} + {V_{ij}\left( {{G_{ij}*{\sin\left( \theta_{ij} \right)}} - {B_{ij}*{\cos\left( \theta_{ij} \right)}}} \right)}}} \end{matrix}$ Note that the variables V_(ij) are given by the Voltage Tensor product represented by bilinear constraints: V _(ij) ={circumflex over (V)} _(i) ·{circumflex over (V)} _(j) In addition, the variables θ_(ij) are given by the Angle Difference relation represented by:

$\theta_{ij} = {\theta_{i} - \theta_{j} - {\sum\limits_{k = 1}^{K_{p}}\;{w_{p\; k}\Delta_{p\; k}}}}$ where

$\sum\limits_{k = 1}^{K_{p}}\;{w_{p\; k}\Delta_{p\; k}}$ is the phase shift due to a Phasor (if any).

Note also that for the sake of convenience for every physical undirected AC branch an OPF model according to an embodiment of the invention model has two logical directed AC branches: (i,j) and (j,i). These branches are not coupled but their parameters are calculated based on the same physical branch characteristics, such as impedance.

The set of Transformer constraints includes:

Joint mixed-integer-linear tap voltage box constraints: w _(tk) ·V _(tk) ≦{circumflex over (V)} _(tk) ≦w _(tk) · V _(tk).

Link with Bus Voltage:

$V_{i} = {\sum\limits_{k = 1}^{K_{t}}\;{{\hat{V}}_{tk} \cdot {t_{tk}.}}}$

Link with Branch (Tapped) Voltage:

${\hat{V}}_{i} = {\sum\limits_{k = 1}^{K_{t}}\;{{\hat{V}}_{tk}.}}$

Note that in a model according to an embodiment of the invention there are two logical transformers corresponding to one physical transformer: one for (I,j) line and the other for (j,i) line. These two logical transformers not only have their parameters generated by the same underlying set of taps of the physical transformer but they also are explicitly coupled through an artificial component referred to as Coupling. Coupling holds any pair of component IDs and couples all of their choice variables one-to-one with each other such as: for (t ₁ ,t ₂)εcoupling: w _(t) ₁ _(k) =w _(t) ₂ _(k). Similarly, Phasors are coupled.

Note that the tap choice variables w_(tk) are restricted to be in [0; 1] box but in addition can be restricted to be binary by a Mixed-Integer integrability constraint: w_(t) ₁ _(k)ε{0; 1} These constraints implies using an MIP solver and is expensive computationally while providing more responsive transformer control.

Impedance Adjustment constraints produce a choice of the admittance matrix (G, B) that depends on transformer/phasor tap position. For each Impedance Adjustment record a choice variable is introduced that is one-to-one (i.e. do not require extra integrability constraints) linked to transformer/phasor tap choice through the following linear constraints: w _(am) ≧w _(tk) ·I( t _(am) ≦t _(tk) ≦ t _(am)) The chosen admittance matrix is then given by the sum:

$\left( {G,B} \right) = {\sum\limits_{k = 1}^{K_{a}}\;{w_{ak} \cdot \left( {G,B} \right)_{ak}}}$

The set of Area constraints involve aggregation of branch power flows into area interface flows:

$P_{f} = {\sum\limits_{e \in {I{(f)}}}\; P_{e}}$ $Q_{f} = {\sum\limits_{e \in {I{(f)}}}\;{Q_{e}.}}$

The set of DC bus constraints comprises the fundamental set of network constraints called the Current Flow Conservation set (CFC):

$I_{i} = {{\sum\limits_{e \in {\delta^{+}{(i)}}}\; I_{e}} - {\sum\limits_{e \in {\delta^{-}{(i)}}}\; I_{e}}}$ Note that I_(i) reflects external current injection from AC/DC converters and thus if no converters are attached to the DC node I then the box constraint for I_(i) is: 0≦I_(i)≦0

The set of DC branch constraints comprises the DC version of the Voltage Law: I _(ij) ·R _(ij) =V _(j) −V _(i) Finally, the set of Converter constraints comprises:

${AC}\text{/}{DC}\mspace{14mu}{Voltage}\mspace{14mu}{relation}\text{:}\mspace{14mu}\begin{matrix} {V_{DC} = {\frac{1}{2}{V_{AC} \cdot \left( {{\cos(\alpha)} + {\cos(\mu)}} \right)}}} \\ {V_{DC} = {{V_{AC} \cdot {\cos\left( {\alpha + \mu} \right)}} - {I_{DC} \cdot R_{com}}}} \end{matrix}$ ${Power}\mspace{14mu}{factor}\mspace{14mu}{relation}\text{:}\mspace{14mu}\begin{matrix} {{\cos(R)} = {\frac{1}{2}\left( {{\cos(\alpha)} + {\cos(\mu)}} \right)}} \\ {Q = {P \cdot {\tan(R)}}} \end{matrix}$

-   -   Active power injection from AC side: P=V_(DC)·I_(DC).

Note that the sign of active power injection is opposite for Rectifier and Inverter.

Finally, an OPF model can handle an objective function of arbitrary form. A basecase OPF objective function according to an embodiment of the invention is a generic convex quadratic cost of all buses' generation and load:

${\min{\sum\limits_{i \in {V{({AC})}}}\;\left\lbrack {{C_{ig}^{a}\left( {\hat{P}}_{ig} \right)}^{2} + {C_{il}^{a}\left( {\hat{P}}_{il} \right)}^{2} + {C_{ig}^{r}\left( {\hat{Q}}_{ig} \right)}^{2} + {C_{il}^{r}\left( {\hat{Q}}_{il} \right)}^{2} + {B_{ig}^{a}{\hat{P}}_{ig}} + {B_{il}^{a}{\hat{P}}_{il}} + {B_{ig}^{r}{\hat{Q}}_{ig}} + {B_{il}^{r}{\hat{Q}}_{il}}} \right\rbrack}},$ where each active/reactive load or generation term reflects a linear transformation, that is, a centering and scaling with respect to a given reference value, of the corresponding model variables. Alternatively, another more complex objective function is one that minimizes the active power losses, defined as:

$\min{\sum\limits_{i = 1}^{n}\;{\sum\limits_{j = 1}^{n}\;{G_{ij}\left( {V_{i}^{2} + \left( \frac{V_{j}}{{ri}_{j}} \right)^{2} - {2\frac{V_{i}V_{j}}{r_{ij}}{\cos\left( {\theta_{i} - \theta_{j} + \varphi_{ij}} \right)}}} \right)}}}$ where V_(i), V_(j) and θ_(i), θ_(j) represent the voltages and angles of the i-th and j-th buses, respectively, and G_(ij) represents the conductances of the branch connecting buses i and j. The ratio of the branch (i,j) is denoted by r_(ij), and the additional angle difference between buses i and j is denoted by φ_(ij). Note that φ_(ij)=0 in the case where there is no phase shifter transformer in the branch (i, j). According an embodiment of the invention, a quadratic approximation of a non-linear objection function will be used.

A basecase OPF model according to an embodiment of the invention as defined above is very complex for numerical optimization due to:

a) Non-convex Voltage Law constraints;

b) Fully coupled Power Flow Conservation network constraints;

c) Integer choice variables for transformer/phasor taps & switching shunts;

d) Large network size in realistic examples;

From the above presentation it is evident that solving OPF models is a very challenging task. Its solution demands the development of sophisticated optimization algorithms which can handle both the large number of decision variables and the complex (non-convex) type of constraints. However, embodiments of the invention provide algorithms that have been successfully applied to the solution of very large scale electrical networks and the smart grid.

SUMMARY

Exemplary embodiments of the invention as described herein generally include systems and methods for approximating an optimal power flow of a smart electric power grid.

According to an aspect of the invention, there is provided a method for approximating an optimal power flow of a smart electric power grid, including providing a cost function that models a smart electric power grid having buses iε[1,N] connected by branches ijε[1,N]×[1,N] with connection matrix

$R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ in terms of bus control vector u_(i) that include voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) that include net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij), that includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij), that includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, constraining the bus and branch control vectors, and the bus and branch dependent variables vectors based on physical and operational requirements and actual measurements, deriving a set of linear equations that minimize the cost function subject to the constraints, from an expression of an extremum of the cost function with respect to all arguments, reducing a dimension of the linear equations by solving for a subset of the linear equations, re-organizing the reduced dimension linear equations into primal and dual parts, and decomposing the re-organized reduced dimensional linear equations into two systems of block matrix equations, where a solution of the two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through the smart electric power grid.

According to a further aspect of the invention, the cost function is

${{f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\;\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}u_{i}^{T}V_{i}^{U}u_{i}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}u_{i}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}\;{R_{ij} \cdot \left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}u_{ij}^{T}W_{ij}^{U}u_{ij}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}u_{ij}}} \right\rbrack}}}},$ where bus control vector u_(i) includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, and connection matrix

$R_{ij} = \left\{ {{{\begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix}\mspace{14mu}{where}\mspace{14mu} x_{i}} = {\sum\limits_{j = 1}^{N}\;{R_{ij} \cdot x_{ij}}}},} \right.$ where V^(X) is a Hessian of the objective function with respect to x_(i), V^(U) is a Hessian of the objective function with respect to u_(i), S^(X) is a gradient of the objective function with respect to x_(i), S^(U) is a gradient of the objective function with respect to u_(i), W^(X) is a Hessian of the objective function with respect to x_(ij), W^(U) is a Hessian of the objective function with respect to u_(ij), T^(X) is a gradient of the objective function with respect to x_(ij), and T^(U) is a gradient of the objective function with respect to u_(ij), subject to the constraints x _(i) ^(LB) ≦x _(i) ≦x _(i) ^(UB) , x _(ij) ^(LB) ≦x _(ij) ≦x _(ij) ^(UB), and u _(i) ^(LB) ≦u _(i) ≦u _(i) ^(UB) , u _(ij) ^(LB) ≦u _(ij) ≦u _(ij) ^(UB), satisfying a voltage law x_(ij)=R_(ij)·(A_(ij) ⁰u_(i)+B_(ij) ⁰u_(j)+C_(ij) ⁰u_(ij)+D_(ij) ⁰u_(ji)+E_(ij) ⁰), where the matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰, u_(j) ⁰, u_(ij) ⁰).

According to a further aspect of the invention, the linear equations that minimize the cost function are:

${{\begin{bmatrix} H_{11} & 0 & {- I} & H_{41}^{T} & H_{51}^{T} & {- I} & 0 \\ 0 & 0 & {- I} & 0 & 0 & 0 & {- I} \\ {- I} & {- I} & 0 & 0 & 0 & 0 & 0 \\ H_{41} & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{51} & 0 & 0 & 0 & 0 & 0 & 0 \\ {{diag}({Zx})} & 0 & 0 & 0 & 0 & {{diag}(X)} & 0 \\ 0 & {{diag}({Zs})} & 0 & 0 & 0 & 0 & {{diag}(S)} \end{bmatrix} \cdot \begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; Z_{x}} \\ {\Delta\; Z_{s}} \end{bmatrix}} = {- \begin{bmatrix} {{\nabla{f(X)}} - Z_{x} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \\ {{Z_{x} \cdot X} - {\mu\; e}} \\ {{Z_{e} \cdot S} - {\mu\; e}} \end{bmatrix}}},$ where [ΔX, ΔS, ΔY, Δy₁, Δy₂, ΔZ_(x), ΔZ_(s)] are unknowns,

$\mspace{20mu}{{H_{11} = {{diag}\left( {V^{x},W^{x},V^{u},W^{U}} \right)}},\mspace{20mu}{H_{41} = \left\lbrack {{- I},R,0,0} \right\rbrack},{R = \begin{pmatrix} \left\lbrack {R_{11}I\mspace{14mu}\ldots\mspace{14mu} R_{1\; N_{1}}I} \right\rbrack & 0 & \ldots & 0 \\ 0 & \left\lbrack {R_{21}I\mspace{14mu}\ldots\mspace{14mu} R_{2\; N_{2}}I} \right\rbrack & \ldots & 0 \\ 0 & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & \left\lbrack {R_{N_{1}1}I\mspace{14mu}\ldots\mspace{14mu} R_{N_{1}N_{1}}I} \right\rbrack \end{pmatrix}},}$ a row of branch (i, j) of H₅₁ is

${\left( H_{51} \right)_{ij} = \left\lbrack {0,{- I},{\ldots\mspace{14mu} R_{ij}A_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}B_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}C_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}D_{ij}^{0}},{\ldots\mspace{14mu} 0}} \right\rbrack},{X = \left( {\left\{ x_{il} \right\},\left\{ x_{ijl} \right\},\left\{ u_{il} \right\},\left\{ u_{ijl} \right\}} \right)^{T}},{{{slack}\mspace{14mu}{variables}\mspace{14mu} S} = \left( {\left\{ s_{xi} \right\},\left\{ s_{xij} \right\},\left\{ s_{ui} \right\},\left\{ s_{uij} \right\}} \right)^{T}},{{o\left( {x_{ij},x_{i}} \right)} = {x_{i\; 1} - {\sum\limits_{j = 1}^{N}\;{R_{ij} \cdot x_{{ij}\; 1}}} - b_{x}}},{{v(X)} = {x_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot u_{{ij}\; 1}} + {D_{ij}^{0} \cdot u_{{ji}\; 1}}} \right)} - b_{v}}},$ Y, y₁, and y₂ are dual variables, and Z_(X)=μ/X, Z_(S)=μ/S where μ is a barrier parameter.

According to a further aspect of the invention, reducing a dimension of the linear equations further comprises solving for ΔZ _(X)=diag(X)⁻¹(−r _(Z) _(X) −diag(Z _(X))ΔX), where ΔZ _(S)=diag(S)⁻¹(−r _(Z) _(S) −diag(Z _(S))ΔS), r _(Z) _(X) =Z _(X) ·X−μe, r _(Z) _(S) =Z _(S) ·S−μe.

According to a further aspect of the invention, the reduced dimension linear equations re-organized into primal and dual parts are

${\begin{bmatrix} H_{41} & 0 & | & 0 & 0 & 0 \\ H_{51} & 0 & | & 0 & 0 & 0 \\ {- I} & {- I} & | & 0 & 0 & 0 \\  - & - & \; & - & - & - \\ {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{x} \right)}}} & 0 & | & {- I} & H_{41}^{T} & H_{51}^{T} \\ 0 & {{{diag}(S)}^{- 1}{{diag}\left( Z_{s} \right)}} & | & {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\  - \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix}} = {\quad{\begin{bmatrix} {- r_{y_{1}}} \\ {- r_{y_{2}}} \\ {- r_{Y}} \\  - \\ {{- r_{X}} - {{{diag}(X)}^{- 1}{r_{Z}.}}} \\ {{- r_{S}} - {{{diag}(S)}^{- 1}{r_{Z}.}}} \end{bmatrix},\mspace{20mu}{{{where}\mspace{20mu}\begin{bmatrix} r_{X} \\ r_{S} \\ r_{Y} \\ r_{y\; 1} \\ r_{y\; 2} \end{bmatrix}} = {\begin{bmatrix} {{\nabla{f(X)}} - Z_{x} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \end{bmatrix}.}}}}$

According to a further aspect of the invention, the two systems of block matrix equations are

${{\begin{bmatrix} H_{41} & 0 \\ H_{51} & 0 \\ {- I} & {- I} \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} = \begin{bmatrix} {- r_{y\; 1}} \\ {- r_{y\; 2}} \\ {- r_{Y}} \end{bmatrix}},{{{{{and}\begin{bmatrix} {\overset{\sim}{H}}_{11} & 0 \\ 0 & {\overset{\sim}{H}}_{22} \end{bmatrix}}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} + {\begin{bmatrix} {- I} & H_{41}^{T} & H_{51}^{T} \\ {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix}}} = \begin{bmatrix} {- {\overset{\sim}{r}}_{X}} \\ {- {\overset{\sim}{r}}_{S}} \end{bmatrix}},{where}$ ${{\overset{\sim}{H}}_{11} = {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{X} \right)}}}},{{\overset{\sim}{H}}_{22} = {{{diag}(S)}^{- 1}{{diag}\left( Z_{S} \right)}}},{{\overset{\sim}{r}}_{X} = {r_{X} + {{{diag}(X)}^{- 1}r_{Z_{X}}}}},{{\overset{\sim}{r}}_{S} = {r_{S} + {{{diag}(S)}^{- 1}r_{Z_{S}}}}},{{\Delta\; X} = \begin{bmatrix} \left\{ {\Delta\; x_{i}} \right\} \\ \left\{ {\Delta\; x_{ij}} \right\} \\ \left\{ {\Delta\; u_{i}} \right\} \\ \left\{ {\Delta\; u_{ij}} \right\} \end{bmatrix}},{{\Delta\; y_{1}} = \left\{ {\Delta\; y_{1}^{x_{i}}} \right\}},{{\Delta\; y_{2}} = {\left\{ {\Delta\; y_{2}^{x_{ij}}} \right\}.}}$

According to a further aspect of the invention, the first block matrix equation is expanded for each bus i and branch(i, j) as follows:

${{\Delta\; x_{i}} - {\sum\limits_{j}\;{R_{ij}\Delta\; x_{ij}}}} = r_{y\; 1}^{x_{i}}$ Δ x_(ij) − R_(ij)(A_(ij)Δ u_(i) + B_(ij)Δ u_(j) + C_(ij)Δ u_(ij) + D_(ij)Δ u_(ji)) = r_(y 2)^(x_(ij)) Δ X = Δ S = rY, and the second block matrix equation is expanded for each bus i and branch(i, j) as follows:

${{{\overset{\sim}{V}}_{i}^{X}\Delta\; x_{i}} - {\Delta\; Y^{x_{i}}} - {\Delta\; y_{1}^{x_{i}}}} = {- {\overset{\sim}{r}}_{x_{i}}}$ ${{{\overset{\sim}{W}}_{ij}^{X}\Delta\; x_{ij}} - {\Delta\; Y^{x_{ij}}} + {R_{ij}\Delta\; y_{1}^{x_{i}}} - {\Delta\; y_{2}^{x_{ij}}}} = {- {\overset{\sim}{r}}_{x_{ij}}}$ ${{{\overset{\sim}{V}}_{i}^{U}\Delta\; u_{i}} - {\Delta\; Y^{u_{i}}} + {\sum\limits_{j}{R_{ij}A_{ij}^{T}\Delta\; y_{2}^{x_{ij}}}} + {\sum\limits_{j}{R_{ji}B_{ji}^{T}\Delta\; y_{2}^{x_{ji}}}}} = {- {\overset{\sim}{r}}_{u_{i}}}$ ${{{\overset{\sim}{W}}_{ij}^{U}\Delta\; u_{ij}} - {\Delta\; Y^{u_{ij}}} + {R_{ij}C_{ij}^{T}\Delta\; y_{2}^{x_{ij}}} + {R_{ji}D_{ji}^{T}\Delta\; y_{2}^{x_{ji}}}} = {- {\overset{\sim}{r}}_{u_{ij}}}$ ${{{{\overset{\sim}{H}}_{22}\Delta\; S} - {\Delta\; Y}} = {- {\overset{\sim}{r}}_{S}}},$ where the method further comprises solving for the other variables in the block matrix equations in the following order: {Δu _(i) },{Δu _(ij) }→{Δx _(ij) }→{Δx _(i) }→S→ΔY→Δy ₁ →Δy ₂.

According to a further aspect of the invention, the method includes initializing unknowns [ΔX^(k), ΔS^(k), ΔY^(k), Δy₁ ^(k), Δy₂ ^(k)] where k is an iteration counter, updating the barrier parameter as μ^(k+1)←ε·μ^(k), εε(0.1, . . . , 0.9), evaluating (ΔX^(k), ΔS^(k), ΔY^(k), Δy_(i) ^(k), Δy₂ ^(k)) using the two block matrix equations, calculating (ΔZ_(X) ^(k), ΔZ_(S) ^(k)), determining a step length Δβ^(k) from a line search, and updating X ^(k+1) ←X ^(k)+Δβ^(k) ΔX ^(k), S ^(k+1) ←S ^(k)+Δβ^(k) ΔS ^(k), Y ^(k+1) ←Y ^(k)+Δβ^(k) ΔY ^(k), y ₁ ^(k+1) ←y ₁ ^(k)+Δβ^(k) Δy ₁ ^(k), Z _(X) ^(k+1) ←Z _(X) ^(k)+Δβ^(k) ΔZ _(X) ^(k), Z _(S) ^(k+1) ←Z _(S) ^(k)+Δβ^(k) ΔZ _(S) ^(k), where the steps of updating the barrier parameter, evaluating (ΔX^(k), ΔS⁴, ΔY^(k), Δy_(i) ^(k), Δy₂ ^(k)), calculating (ΔZ_(X) ^(k), ΔZ_(S) ^(k)), determining a step length Δβ^(k), and updating X^(k+1), S^(k+1), Y^(k+1), y₁ ^(k+1), y₂ ^(k+1), Z_(X) ^(k+1) and Z_(S) ^(k+1) are repeated until either a dual gap or a maximum of (r_(X), r_(S)) is less than a pre-specified minimum.

According to a further aspect of the invention, the cost function is,

${{f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}{\sum\limits_{Ki}{u_{i}^{k\; T}V_{i}^{U}u_{i}^{k}}}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}{\sum\limits_{Ki}u_{i}^{k}}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}{R_{ij} \cdot \left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}{\sum\limits_{Kij}{u_{ij}^{k\; T}W_{ij}^{U}u_{ij}^{k}}}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}{\sum\limits_{Kik}u_{ij}^{k}}}} \right\rbrack}}}},$ where bus control vector u_(i) includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, and connection matrix

$R_{ij} = \left\{ {{{\begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix}{where}\mspace{14mu} x_{i}} = {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{ij}}}},} \right.$ where V^(X) is a Hessian of the objective function with respect to x_(i), V^(U) is a Hessian of the objective function with respect to u_(i), S^(X) is a gradient of the objective function with respect to x_(i), S^(U) is a gradient of the objective function with respect to u_(i) W^(X) is a Hessian of the objective function with respect to x_(ij), W^(U) is a Hessian of the objective function with respect to u_(ij), T^(X) is a gradient of the objective function with respect to x_(ij), and T^(U) is a gradient of the objective function with respect to u_(ij), subject to constraints

${x_{i}^{LB} \leq x_{i} \leq x_{i}^{UB}},{u_{i}^{LB} \leq {\sum\limits_{Ki}u_{i}^{k}} \leq u_{i}^{UB}},{x_{ij}^{LB} \leq x_{ij} \leq x_{ij}^{UB}},{u_{ij}^{LB} \leq {\sum\limits_{Kij}u_{ij}^{k}} \leq u_{ij}^{UB}},$ where a control grid is discretized by binary variables zε{0,1} as z _(i) ^(k) u _(ik) ^(LB) ≦u _(i) ^(k) ≦z _(i) ^(k) u _(ik) ^(UB), z _(ij) ^(k) u _(ijk) ^(LB) ≦u _(ij) ^(k) ≦z _(ij) ^(k) u _(ijk) ^(UB), where the binary variables are subject to constraints

$z_{i}^{LB} \leq {\sum\limits_{Ki}z_{i}^{k}} \leq z_{i}^{UB}$ $z_{i}^{LB} \leq {\sum\limits_{Ki}z_{i}^{k}} \leq z_{i}^{UB}$ and satisfying a voltage law x_(ij)=R_(ij)·(Σ_(Ki)A_(ij) ^(k)u_(i) ^(k)+Σ_(Kj)B_(ij) ^(k)u_(j) ^(k)+Σ_(Kij)C_(ij) ^(k)u_(ij) ^(k)+Σ_(Kji)D_(ij) ^(k)u_(ji) ^(k)) where the matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰, u_(j) ⁰, u_(ij) ⁰).

According to a further aspect of the invention, the cost function is ƒ({tilde over (x)},u,y,v)=½[{tilde over (x)} ^(T) V ^(X) {tilde over (x)}+u ^(T) V ^(U) u+y ^(T) W ^(Y) y+v ^(T) W ^(V) v]+S ^(X) {tilde over (x)}+S ^(U) u+T ^(Y) y+T ^(V) v, where bus control vector u includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x includes net active X_(i) ^(P) and reactive X_(i) ^(Q) power injection, branch control vector v includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector y includes net active X_(ij) ^(P) and reactive X_(ij) ^(Q) power injections, where bus dependent vector is defined as

${x = {\left. {\overset{\sim}{A}\overset{\sim}{x}}\Leftrightarrow x_{i} \right. = {{\sum\limits_{k = 1}^{G_{i}}x_{i,k}^{g}} - {\sum\limits_{k = 1}^{L_{i}}x_{i,k}^{l}} + {\sum\limits_{k = 1}^{{SH}_{i}}x_{i,k}^{sh}}}}},$ where G_(i), L, and SH_(i) are the number of generators, loads and shunts placed on bus i and x_(i,k) ^(g), x_(i,k) ^(l) and x_(i,k) ^(sh) are the k^(th) active and reactive generation, load and shunt variables at bus i, respectively, and Ãε{−1,0,1}^(N×2(G+L+SH)), where V^(X) is a Hessian of the objective function with respect to x, V^(U) is a Hessian of the objective function with respect to u, S^(X) is a gradient of the objective function with respect to x, S^(U) is a gradient of the objective function with respect to u, W^(Y) is a Hessian of the objective function with respect to y, W^(V) is a Hessian of the objective function with respect to v, T^(Y) is a gradient of the objective function with respect to y, and T^(V) is a gradient of the objective function with respect to v, the cost function having connection matrix

$R_{ij} = \left\{ {{\begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix}{where}\mspace{14mu} a_{i}\overset{\sim}{x}} = {\sum\limits_{j = 1}^{N}{R_{ij}y_{ij}}}} \right.$ where a_(i) is i^(th) is row of matrix Ã, subject to the constraints {tilde over (x)} ^(LB) ≦{tilde over (x)}≦{tilde over (x)} ^(UB) , y ^(LB) ≦y≦y ^(UB), u ^(LB) ≦u≦u ^(UB) , v ^(LB) ≦v≦v ^(UB), and satisfying a voltage law y_(ij)=R_(ij)·(A_(ij) ⁰·u_(i)+B_(ij) ⁰·u_(j)+C_(ij) ⁰·v_(ij)+D_(ij) ⁰·v_(ji)+E_(ij) ⁰), where the matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰, u_(j) ⁰, u_(ij) ⁰).

According to another aspect of the invention, there is provided a method for approximating an optimal power flow of a smart electric power grid, including providing a cost function that models a smart electric power grid having buses iε[1,N] connected by branches ijε[1, N]×[1, N] with connection matrix

$R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ in terms of bus control vector u_(i) that include voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) that include net active X_(i) ^(P) and reactive X_(i) ^(Q) power injection, branch control vector u_(ij), that includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij), that includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, constraining the bus and branch control vectors, and the bus and branch dependent variables vectors based on physical and operational requirements and actual measurements, deriving a set of non-linear vector equations

$\overset{\sim}{L} = {\begin{bmatrix} {{\nabla{f(X)}} - {\mu/X} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{{- \mu}/S} - Y} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \end{bmatrix} = 0}$ that minimize the cost function subject to the constraints, from an expression of an extremum of the cost function with respect to all arguments, where X=({x_(i)}, {x_(ij)}, {u_(i)}, {u_(ij)})^(T), bus control vector u_(i), includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x, includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, S=({s_(xi)}, {s_(xij)}, {s_(ui)}, {s_(uij)})^(T), where s_(xi), s_(ui), s_(xij), s_(uij) are slack variables,

${{o\left( {x_{ij},x_{i}} \right)} = {x_{i\; 1} - {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{{ij}\; 1}}} - b_{x}}},{{v(X)} = {x_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot u_{{ij}\; 1}} + {D_{ij}^{0} \cdot u_{{ji}\; 1}}} \right)} - b_{v}}}$ b_(x) and b_(v) are predetermined constants,

$R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ is a connection matrix where

${x_{i\;} = {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{{ij}\;}}}},{x_{{ij}\;} = {R_{ij} \cdot \left( {{A_{ij}^{0}u_{i\;}} + {B_{ij}^{0}u_{j\;}} + {C_{ij}^{0}u_{{ij}\;}} + {D_{ij}^{0}u_{{ji}\;}} + E_{ij}^{0}} \right)}}$ is a voltage law where the matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰, u_(j) ⁰, u_(ij) ⁰), Y, y₁, and y₂ are dual variables, and μ is a barrier parameter, re-organizing the vector {tilde over (L)} as

${\overset{\sim}{L} = \begin{bmatrix} \left\{ L^{i} \right\} \\ \left\{ L^{ij} \right\} \end{bmatrix}},$ where each L^(i) represent the constraints of bus i, and each L^(ij) represents the constraints of branch (i, j), applying Newton's method to every bus L^(i) and branch L^(ij), to obtain

${{{H^{i} \cdot \begin{bmatrix} {\Delta\; X^{i}} \\ {\Delta\; S^{i}} \\ {\Delta\; Y^{i}} \\ {\Delta\; y_{1}^{i}} \\ {\Delta\; y_{2}^{i}} \end{bmatrix}} = {- L^{i}}};{{H^{ij} \cdot \begin{bmatrix} {\Delta\; X^{ij}} \\ {\Delta\; S^{ij}} \\ {\Delta\; Y^{ij}} \\ {\Delta\; y_{1}^{ij}} \\ {\Delta\; y_{2}^{ij}} \end{bmatrix}} = {- L^{ij}}}},$ where H^(i) and H^(ij) are Jacobian matrices, ΔX^(i)=({Δx_(i)}, {Δu_(i)})^(T), ΔX^(ij)=({Δx_(ij)}, {Δu_(ij)}), ΔS^(i)=({Δs_(xi)}, {Δs_(ui)})^(T), ΔS^(ij)=({Δs_(xij)}, {Δs_(uij)}), and [ΔX, ΔS, ΔY, Δy₁, Δy₂] are unknowns, and using cutting planes to solve these block-wise Newton's equations, where a solution of the two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through the smart electric power grid.

According to another aspect of the invention, there is provided a non-transitory program storage device readable by a computer, tangibly embodying a program of instructions executed by the computer to perform the method steps for approximating an optimal power flow of a smart electric power grid.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram of an exemplary, non-limiting electric power system network, according to an embodiment of the invention.

FIG. 2 is a flow chart of an interior point method for solving an optimal power flow model of a smart power grid, according to an embodiment of the invention.

FIG. 3 is a flow chart of an iterative method for solving an optimal power flow model of a smart power grid, according to an embodiment of the invention.

FIG. 4 is a block diagram of an exemplary computer system for implementing an interior point method for solving an optimal power flow model of a smart power grid, according to an embodiment of the invention.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

Exemplary embodiments of the invention as described herein generally include systems and methods for estimating a solution of an optimal power flow model of a smart power grid. Accordingly, while the invention is susceptible to various modifications and alternative forms, specific embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that there is no intent to limit the invention to the particular forms disclosed, but on the contrary, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention.

An optimal power flow (OPF) model seeks to produce an optimal power flow solution that balances generation/consumption for all buses with power distribution across all connected branches such that the overall system state, as measured by underlying physical parameters including voltages and phase angles, is feasible and safe. A basecase OPF model can be used to control AC and DC energy dispatch (i.e. bus generation/load and branch power flow), regulate AC bus voltages and phase angles, select AC transformer/phasor taps and shunt switches, and control an AC/DC converter.

Solver Model Specification

An electric grid Optimal Power Flow (OPF) model according to an embodiment of the invention includes buses iε[1,N] connected by branches ijε[1,N]×[1,N]. The AC-AC branch (i.e., a branch obeying alternating current physics laws) connection matrix R is given as:

$R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {o/{w.}} \end{matrix} \right.$ Minimal AC OPF Model Specification

A minimal OPF model according to an embodiment of the invention assumes an AC-only grid with continuous-only state-dependent variables x and controls u. In addition to that, there are no thermal capacity limits enforced on branch power flows. Finally, discrete logic z defining a linearization region of the originally non-convex AC voltage law constraints is assumed to be a fixed binary vector z₀.

Box-Constraints:

According to an embodiment of the invention, the following two box-constraint families are defined that independently limit bus control vector u_(i) and dependent bus variable vector x_(i): x _(i) ^(LB) ≦x _(i) ≦x _(i) ^(UB), u _(i) ^(LB) ≦u _(i) ≦u _(i) ^(UB),  (1) where the superscripts LB and UB refer to lower bound and upper bound, respectively. Bus control vector u_(i) typically includes voltage magnitude V_(i), log-scaled in this model, phase angle θ_(i), and shunt switch position Q_(i•). Bus dependent variables vector x_(i) typically includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection. Similarly, according to an embodiment of the invention, a box-constraint family is defined for branch variables that independently limit branch control vector u_(ij), which typically includes transformer tap selection t_(ij) and phase shift φ_(ij), and dependent branch variable vector x_(ij), which typically includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections: x _(ij) ^(LB) ≦x _(ij) ≦x _(ij) ^(UB), u _(ij) ^(LB) ≦u _(ij) ≦u _(ij) ^(UB),  (2)

Flow Balance:

At every bus a net power injection balance should be maintained between a bus and connected branches. Branch dependent variables vector x_(ij), in which typically x_(ij)=(x_(ij) ^(P), x_(ij) ^(Q)), can be defined to specify this balance as follows:

$\begin{matrix} {x_{i} = {\sum\limits_{j = 1}^{N}\;{R_{ij} \cdot {x_{ij}.}}}} & (3) \end{matrix}$

The matrix R is always highly sparse since the connection density of real-world power grids is typically 1-5 branches per bus.

Discrete Control Grid:

Some of the components in bus control vector u_(i) may be discrete, while others may be continuous. According to an embodiment of the invention, there is a minimal discretization of the overall Cartesian product space of the vector u_(i) into intervals kε[1,K_(i)] each defined by a box [u_(ik) ^(LB), u_(ik) ^(LB)], a fixed anchor vector u_(ik) ⁰ and a binary choice variable z_(ik). The latter is assumed fixed at z_(ik) ⁰ in this minimal AC OPF specification. Then, a discrete control grid is given by the following discrete box-constraint family:

$\begin{matrix} {{\sum\limits_{k = 1}^{K_{i}}\;{Z_{ik}^{0}u_{ik}^{LB}}} \leq u_{i} \leq {\sum\limits_{k = 1}^{K_{i}}\;{Z_{ik}^{0}{u_{ik}^{UB}.}}}} & (4) \end{matrix}$ Similarly, according to an embodiment of the invention, a discrete control grid is defined for the branch control vector u_(ij) defined by kε└1,K_(ij)┘ intervals each defined by a box [u_(ijk) ^(LB), u_(ijk) ^(UB)], a fixed anchor u_(ijk) ⁰ and a binary choice variable z_(ijk), assumed fixed at z_(ijk) ⁰ in this minimal AC OPF specification:

$\begin{matrix} {{\sum\limits_{k = 1}^{K_{ji}}\;{Z_{ijk}^{0}u_{ijk}^{LB}}} \leq u_{ij} \leq {\sum\limits_{k = 1}^{K_{ij}}\;{Z_{ijk}^{0}{u_{ijk}^{UB}.}}}} & (5) \end{matrix}$

Discretized Voltage Law:

Each interval k implies a new feasibility set X_(ij) ^(k) for the branch power flow vector x_(ij) represented as a rectangle generated by intervals [u_(ik) ^(LB), u_(ik) ^(UB)], [u_(jk) ^(LB), u_(jk) ^(UB)], and [u_(ijk) ^(LB), u_(ijk) ^(UB)], following the voltage law gradient lines drawn at the combined anchor point (u_(ik) ⁰, u_(jk) ⁰, u_(ijk) ⁰). According to an embodiment of the invention, the feasibility set X_(ij) ^(k) is not explicitly included in the model specification. Instead, a discretized version of the AC voltage law is provided that defines the relation between branch power flow vectors x_(ij) and control vectors u_(i), u_(j), u_(ij) that are linearized with respect to a given anchor point (u_(i) ⁰, u_(j) ⁰, u_(ij) ⁰) that is typically given by the grid anchor point (u_(ik) ⁰, u_(jk) ⁰, u_(ijk) ⁰) corresponding to z_(ik) ⁰ and z_(ijk) ⁰. Then, an embodiment of the invention assumes the following form of the discretized AC voltage law: x _(ij) =R _(ij)·(A _(ij) ⁰ u _(i) +B _(ij) ⁰ u _(j) +C _(ij) ⁰ u _(ij) +D _(ij) ⁰ u _(ji) +E _(ij) ⁰).  (6) Note that this constraint exhibits two types of sparsity: one inherited from R_(ij) and the other due to dependence on only two control vectors (u_(i), u_(j)) and (u_(ij), u_(ji)) instead of the whole-grid control vector U=└(u_(i))_(i=1) ^(N), (u_(ij))_(i,j=1) ^(N)┘. The matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for the given anchor point (u_(i) ⁰, u_(j) ⁰, u_(ij) ⁰).

Objective Function:

According to an embodiment of the invention, it is assumed that there exists a generic convex quadratic cost function to be minimized. This function includes terms related to managing the control vector u_(i) that minimizes control adjustment costs, bus dependent variables vector x_(i) that minimize generation costs, and branch dependent variables vector x_(ij) that minimize transmission losses:

$\begin{matrix} {{{f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\;\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}u_{i}^{T}V_{i}^{U}u_{i}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}u_{i}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}\;{R_{ij} \cdot \left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}u_{ij}^{T}W_{ij}^{U}u_{ij}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}u_{ij}}} \right\rbrack}}}},} & (7) \end{matrix}$ where:

V^(X) is a Hessian of the objective function with respect to x_(i),

V^(U) is a Hessian of the objective function with respect to u_(i),

S^(X) is a gradient of the objective function with respect to x_(i),

S^(U) is a gradient of the objective function with respect to u_(i),

W^(X) is a Hessian of the objective function with respect to x_(ij),

W^(U) is a Hessian of the objective function with respect to u_(ij),

T^(X) is a gradient of the objective function with respect to x_(ij), and

T^(U) is a gradient of the objective function with respect to u_(ij).

Discrete AC OPF Model Specification

A discrete OPF model according to an embodiment of the invention assumes that a discrete logic z defines a linearization region of the originally non-convex AC voltage law constraints is also a model variable that is constrained to be a binary zε{0,1}. A discrete OPF model according to an embodiment of the invention has a form similar to the minimal OPF except for the following changes:

Box-Constraints:

The dependent variables box-bounds remain unchanged: x _(i) ^(LB) ≦x _(i) ≦x _(i) ^(UB), x _(ij) ^(LB) ≦x _(ij) ≦x _(ij) ^(UB).  (8) However, the box-bound constraints for control variables u_(i) and u_(ij) are changed to have the following simplex form:

$\begin{matrix} {{u_{i}^{LB} \leq {\sum\limits_{Ki}\; u_{i}^{k}} \leq u_{i}^{UB}}{u_{ij}^{LB} \leq {\sum\limits_{Kij}\; u_{ij}^{k}} \leq {u_{ij}^{UB}.}}} & (9) \end{matrix}$ Also, a simplex form is introduced for the discrete variables:

$\begin{matrix} {{z_{i}^{LB} \leq {\sum\limits_{Ki}\; z_{i}^{k}} \leq z_{i}^{UB}}{z_{i}^{LB} \leq {\sum\limits_{Ki}\; z_{i}^{k}} \leq z_{i}^{UB}}} & (10) \end{matrix}$ Discrete Control Grid: z _(i) ^(k) u _(ik) ^(LB) ≦u _(i) ^(k) ≦z _(i) ^(k) u _(ik) ^(UB), z _(ij) ^(k) u _(ijk) ^(LB) ≦u _(ij) ^(k) ≦z _(ij) ^(k) u _(ijk) ^(UB).  (11)

Note that the vector u_(ij) ^(k) of branch control variables, constrained within a discrete control grid interval k, now also includes power flow anchor points (P_(ij) ^(k), Q_(ij) ^(k)) for every branch.

Discretized Voltage Law. x _(ij) =R _(ij)·(Σ_(Ki) A _(ij) ^(k) u _(i) ^(k)+Σ_(Ki) B _(ij) ^(k) u _(j) ^(k)+Σ_(Kij) C _(ij) ^(k) u _(ij) ^(k)+Σ_(Kji) D _(ij) ^(k) u _(ji) ^(k))  (12)

Objective Function:

$\begin{matrix} {{f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\;\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}{\sum\limits_{Ki}\;{u_{i}^{kT}V_{i}^{U}u_{i}^{k}}}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}{\sum\limits_{Ki}\; u_{i}^{k}}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}\;{R_{ij} \cdot \left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}{\sum\limits_{Kij}\;{u_{ij}^{kT}W_{ij}^{U}u_{ij}^{k}}}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}{\sum\limits_{Kik}\; u_{ij}^{k}}}} \right\rbrack}}}} & (13) \end{matrix}$ Interior Point Method Solver Math for a Minimal AC OPF Model

First, a model according to an embodiment of the invention is standardized to the following format:

$\begin{matrix} {{\min\;{f\left( {x,u} \right)}} = {{\sum\limits_{i = 1}^{N}\;\left\lbrack {{\frac{1}{2}{x_{i\; 1}^{T} \cdot V_{i}^{X} \cdot x_{i\; 1}^{T}}} + {\frac{1}{2}{u_{i\; 1}^{T} \cdot V_{i}^{U}}u_{i\; 1}^{T}} + {S_{i\; 1}^{X} \cdot x_{i\; 1}} + {S_{i\; 1}^{U} \cdot u_{i\; 1}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}\;{R_{ij} \cdot \begin{bmatrix} {{\frac{1}{2}{x_{{ij}\; 1}^{T} \cdot W_{{ij}\; 1}^{X} \cdot x_{{ij}\; 1}^{T}}} +} \\ {{\frac{1}{2}{u_{ij}^{T} \cdot W_{ij}^{U} \cdot u_{{ij}\; 1}^{T}}} + {T_{{ij}\; 1}^{X} \cdot x_{{ij}\; 1}} + {T_{{ij}\; 1}^{U} \cdot u_{{ij}\; 1}}} \end{bmatrix}}}}} & (14.1) \\ {\mspace{79mu}{{X + S} = {X^{UB}\left\{ \begin{matrix} {{x_{i\; 1} + s_{xi}} = x_{i\; 1}^{UB}} \\ {{u_{i\; 1} + s_{ui}} = u_{i\; 1}^{UB}} \\ {{x_{{ij}\; 1} + s_{xij}} = x_{{ij}\; 1}^{UB}} \\ {{u_{{ij}\; 1} + s_{uij}} = u_{{ij}\; 1}^{UB}} \end{matrix} \right.}}} & (14.2) \\ {\mspace{79mu}{{x_{i\; 1} - {\sum\limits_{j = 1}^{N}\;{R_{ij} \cdot x_{{ij}\; 1}}}} = b_{x}}} & (14.3) \\ {\mspace{79mu}{{x_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot u_{{ij}\; 1}} + {D_{ij}^{0} \cdot u_{{ji}\; 1}}} \right)}} = b_{v}}} & (14.4) \\ {\mspace{79mu}{x_{i\; 1},x_{{ij}\; 1},u_{i\; 1},u_{{ij}\; 1},s_{xi},s_{ui},s_{xij},{s_{uij} \geq 0},}} & (14.5) \end{matrix}$ where X=({x_(i1)}, {x_(ij1)}, {u_(u1)}, {u_(ij1)})^(T), s_(xi), s_(ui), s_(xij), s_(uij) are slack variables, S=({s_(xi)}, {s_(xij)}, {s_(ui)}, {s_(uij)})^(T), and

$\begin{matrix} {\mspace{79mu}{{x_{i\; 1} = {x_{i} - x_{i}^{LB}}},\mspace{20mu}{x_{{ij}\; 1} = {x_{ij} - x_{ij}^{LB}}},\mspace{20mu}{x_{i\; 1}^{UB} = {x_{i}^{UB} - x_{i}^{LB}}},\mspace{20mu}{x_{{ij}\; 1}^{UB} = {x_{ij}^{UB} - x_{ij}^{LB}}},}} & (15.1) \\ {\mspace{79mu}{{u_{i\; 1} = {u_{i} - {\max\left\{ {u_{i}^{LB},{\sum\limits_{k = 1}^{K_{i}}\;{z_{ik}^{0} \cdot u_{ik}^{LB}}}} \right\}}}},\mspace{20mu}{u_{{ij}\; 1} = {u_{ij} - {\max\left\{ {u_{ij}^{LB},{\sum\limits_{k = 1}^{K_{ij}}\;{z_{ijk}^{0} \cdot u_{ijk}^{LB}}}} \right\}}}},}} & (15.2) \\ {\mspace{85mu}{{u_{i\; 1}^{UB} = {{\min\left\{ {u_{i}^{UB},{\sum\limits_{k = 1}^{K_{i}}\;{z_{ik}^{0} \cdot u_{ik}^{UB}}}} \right\}} - {\max\left\{ {u_{i}^{LB},{\sum\limits_{k = 1}^{K_{i}}\;{z_{ik}^{0} \cdot u_{ik}^{LB}}}} \right\}}}},\mspace{20mu}{u_{{ij}\; 1}^{UB} = {{\min\left\{ {u_{ij}^{UB},{\sum\limits_{k = 1}^{K_{ij}}\;{z_{ijk}^{0} \cdot u_{ijk}^{UB}}}} \right\}} - {\max\left\{ {u_{ij}^{LB},{\sum\limits_{k = 1}^{K_{ij}}\;{z_{ijk}^{0} \cdot u_{ijk}^{LB}}}} \right\}}}},}} & (15.3) \\ {\mspace{79mu}{{b_{x} = {{\sum\limits_{j = 1}^{N}\;{R_{ij} \cdot x_{ij}^{LB}}} - x_{i}^{LB}}},}} & (15.4) \\ {b_{v} = {{R_{ij} \cdot \left( {{{A_{ij}^{0} \cdot \max}\left\{ {u_{i}^{LB},{\sum\limits_{k = 1}^{K_{i}}\;{z_{ik}^{0} \cdot u_{ik}^{LB}}}} \right\}} + {{B_{ij}^{0} \cdot \max}\left\{ {u_{j}^{LB},{\sum\limits_{k = 1}^{K_{i}}\;{z_{ik}^{0} \cdot u_{jk}^{LB}}}} \right\}} + {{C_{ij}^{0} \cdot \max}\left\{ {u_{ij}^{LB},{\sum\limits_{k = 1}^{K_{ij}}\;{z_{ijk}^{0} \cdot u_{ijk}^{LB}}}} \right\}} + {{D_{ij}^{0} \cdot \max}\left\{ {u_{ij}^{LB},{\sum\limits_{k = 1}^{K_{ij}}\;{z_{ijk}^{0} \cdot u_{jik}^{LB}}}} \right\}} + E_{{ij}\; 1}^{0}} \right)} - {x_{ij}^{LB}.}}} & (15.5) \end{matrix}$ The coefficients of the objective function can also be changed:

$\begin{matrix} {{S_{i\; 1}^{X} = {S_{i}^{X} + {x_{i}^{LB} \cdot V_{i}^{X}}}},{S_{i\; 1}^{U} = {S_{i}^{U} + {\max{\left\{ {u_{j}^{LB},{\sum\limits_{k = 1}^{K_{i}}\;{z_{ik}^{0} \cdot u_{jk}^{LB}}}} \right\} \cdot V_{i}^{U}}}}},{T_{{ij}\; 1}^{X} = {T_{ij}^{X} + {x_{i}^{LB} \cdot W_{ij}^{X}}}},{T_{{ij}\; 1}^{U} = {T_{ij}^{U} + {\max{\left\{ {u_{j}^{LB},{\sum\limits_{k = 1}^{K_{i}}\;{z_{ik}^{0} \cdot u_{jk}^{LB}}}} \right\} \cdot W_{ij}^{U}}}}},} & (16) \end{matrix}$

Let

$\begin{matrix} {{v(X)} = {x_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot u_{{ij}\; 1}} + {D_{ij}^{0} \cdot u_{{ji}\; 1}}} \right)} - {b_{v}.}}} & (17) \\ {{o\left( {x_{ij},x_{i}} \right)} = {x_{i\; 1} - {\sum\limits_{j = 1}^{N}\;{R_{ij} \cdot x_{{ij}\; 1}}} - {b_{x}.}}} & (18) \end{matrix}$

Then, by introducing logarithmic barrier terms, one can obtain the following Lagrangian function: L _(μ)(X,S)=ƒ(X)−μ^(T)(ln X+ln S)−Y ^(T)(X+S−X ^(UB))−y ₁ ^(T) o(x _(ij) ,x _(i))−y ₂ ^(T) v(X))  (19) where X=({x_(i1)}, {x_(ij1)}, {u_(i1)}, {u_(ij1)})^(T), S=({s_(xi)}, {s_(xij)}, {s_(ui)}, {s_(uij)})^(T). The Lagrange multipliers Y and y are called dual variables. According to an embodiment of the invention, the Karush-Kuhn-Tucker (KKT) optimality conditions of the barrier system may be obtained from the gradient of the Lagrangian with respect to all variables, and solving for the extremum:

$\begin{matrix} {\begin{bmatrix} {\nabla L_{X}} \\ {\nabla L_{S}} \\ {\nabla L_{Y}} \\ {\nabla L_{y_{1}}} \\ {\nabla L_{y_{2}}} \end{bmatrix} = {\begin{bmatrix} {{\nabla{f(X)}} - {\mu/X} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{{- \mu}/S} - Y} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \end{bmatrix} = 0.}} & (20) \end{matrix}$ Introducing variables Z_(x)=μ/X, Z_(s)=μ/S transforms EQ. (20) into:

$\begin{matrix} {{\begin{bmatrix} {{\nabla{f\left( {X,Y} \right)}} - Z_{x} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{- Y} - Z_{x}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \\ {{Z_{x} \cdot X} - {\mu\; e}} \\ {{Z_{s} \cdot S} - {\mu\; e}} \end{bmatrix} = 0},} & (21) \end{matrix}$ where e is a unit vector. Newton's method can be used to solve the above non-linear equation, to yield the following linear equation:

$\begin{matrix} {{H \cdot \begin{bmatrix} {\Delta X} \\ {\Delta S} \\ {\Delta Y} \\ {\Delta y}_{1} \\ {\Delta y}_{2} \\ {\Delta Z}_{x} \\ {\Delta Z}_{s} \end{bmatrix}} = {- \begin{bmatrix} {{\nabla{f\left( {X,Y} \right)}} - Z_{x} - Y - {y_{1} \cdot {\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2} \cdot {\nabla{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \\ {{Z_{S} \cdot X} - {\mu\; e}} \\ {{Z_{S} \cdot S} - {\mu\; e}} \end{bmatrix}}} & (22) \end{matrix}$ in which H is the following Jacobian matrix:

$\begin{matrix} {{H = \begin{bmatrix} H_{11} & 0 & {- I} & H_{41}^{T} & H_{51}^{T} & {- I} & 0 \\ 0 & 0 & {- I} & 0 & 0 & 0 & {- I} \\ {- I} & {- I} & 0 & 0 & 0 & 0 & 0 \\ H_{41} & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{51} & 0 & 0 & 0 & 0 & 0 & 0 \\ {{diag}\left( Z_{X} \right)} & 0 & 0 & 0 & 0 & {{diag}(X)} & 0 \\ 0 & {{diag}\left( Z_{S} \right)} & 0 & 0 & 0 & 0 & {{diag}(S)} \end{bmatrix}},} & (23) \end{matrix}$ where H ₁₁=diag(V ^(X) ,W ^(X) v ^(U) ,W ^(U)),  (24.1) where diag( ) is a vector formed from the diagonal components of the matrix argument,

$\begin{matrix} {\mspace{79mu}{{H_{41} = \left\lbrack {{- I},R,0,0} \right\rbrack},}} & (24.2) \\ {R = \begin{pmatrix} \left\lbrack {R_{11}I\mspace{14mu}\ldots\mspace{14mu} R_{{1N_{1}}\;}I} \right\rbrack & 0 & \ldots & 0 \\ 0 & \left\lbrack {R_{21}I\mspace{14mu}\ldots\mspace{14mu} R_{2N_{2}}I} \right\rbrack & \ldots & 0 \\ 0 & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & \left\lbrack {R_{N_{1}1}I\mspace{14mu}\ldots\mspace{14mu} R_{N_{1}N_{1}}I} \right\rbrack \end{pmatrix}} & (24.3) \end{matrix}$ and the row of branch (i, j) of H₅₁ is (H ₅₁)_(ij)=[0,−I, . . . R _(ij) A _(ij) ⁰ , . . . R _(ij) B _(ij) ⁰ , . . . R _(ij) C _(ij) ⁰ , . . . R _(ij) D _(ij) ⁰, . . . 0].  (24.4) Note that, in EQ. (24.1), the arguments V^(X), W^(X), V^(U), W^(U) themselves form a block diagonal

$\begin{matrix} {\mspace{79mu}{{{matrix}{{\text{:}\mspace{14mu}\begin{bmatrix} V^{X} & 0 & 0 & 0 \\ 0 & W^{X} & 0 & 0 \\ 0 & 0 & V^{U} & 0 \\ 0 & 0 & 0 & W^{U} \end{bmatrix}}.\mspace{79mu}{{Define}\mspace{76mu}\begin{bmatrix} r_{X} \\ r_{S} \\ r_{Y} \\ r_{y_{1}} \\ r_{y_{2}} \\ r_{Z_{x}} \\ r_{Z_{s}} \end{bmatrix}}}} = {\begin{bmatrix} {{\nabla{f(X)}} - Z_{x} - Y - {y_{1}^{T}{\nabla\;{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla\;{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- \;{v(X)}} \\ {{Z_{x} \cdot X} - {\mu\; e}} \\ {{Z_{s} \cdot S} - {\mu\; e}} \end{bmatrix}.}}} & (25) \\ {\mspace{76mu}{{Then},}} & \; \\ {{\begin{bmatrix} H_{11} & 0 & {- I} & H_{41}^{T} & H_{51}^{T} & {- I} & 0 \\ 0 & 0 & {- I} & 0 & 0 & 0 & {- I} \\ {- I} & {- I} & 0 & 0 & 0 & 0 & 0 \\ H_{41} & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{51} & 0 & 0 & 0 & 0 & 0 & 0 \\ {{diag}\left( Z_{X} \right)} & 0 & 0 & 0 & 0 & {{diag}(X)} & 0 \\ 0 & {{diag}\left( Z_{S} \right)} & 0 & 0 & 0 & 0 & {{diag}(S)} \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; Z_{x}} \\ {\Delta\; Z_{s}} \end{bmatrix}} = {- \;{\begin{bmatrix} r_{X} \\ r_{S} \\ r_{Y} \\ r_{y_{1}} \\ r_{y_{2}} \\ r_{Z_{x}} \\ r_{Z_{s}} \end{bmatrix}.}}} & (26) \end{matrix}$ Note that the dimension of the system of EQS. (26) can be very large. As a result the solution time needed to solve it may be prohibitive for large scale electrical networks. However, a procedure according to an embodiment of the invention that uses a series of back substitutions, described below, can reduce the dimensionality of the system so that the equations require much less computational time to solve.

According to an embodiment of the invention, to solve these Newton equations, the dimension can be reduced as follows. From the last two equations, one can obtain: ΔZ _(x)=diag(X)⁻¹(−r _(Z) _(x) −diag(Z _(X))ΔX) ΔZ _(s)=diag(S)⁻¹(−r _(Z) _(s) −diag(Z _(S))ΔS)  (27) The solutions of EQS. (27) can be substituted into the other equations of EQS. (26), which can be reorganized into primal and dual parts:

$\begin{matrix} {\begin{bmatrix} H_{41} & 0 & ❘ & 0 & 0 & 0 \\ H_{51} & 0 & ❘ & 0 & 0 & 0 \\ {- I} & {- I} & ❘ & 0 & 0 & 0 \\  - & - & \; & - & - & - \\ {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{x} \right)}}} & 0 & ❘ & {- I} & H_{41}^{T} & H_{51}^{T} \\ 0 & {{{diag}(S)}^{- 1}{{diag}\left( Z_{s} \right)}} & ❘ & {- I} & 0 & 0 \end{bmatrix}{\quad{\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\  - \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix} = {\quad\begin{bmatrix} {- r_{y\; 1}} \\ {- r_{y\; 2}} \\ {- r_{Y}} \\  - \\ {{- r_{x}} - {{{diag}(X)}^{- 1}r_{Z_{x}}}} \\ {{- r_{s}} - {{{diag}(S)}^{- 1}r_{Z_{s}}}} \end{bmatrix}}}}} & (28) \end{matrix}$ By defining

$\begin{matrix} {{{\overset{\sim}{H}}_{11} = {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{x} \right)}}}}{{\overset{\sim}{H}}_{22} = {{{diag}(S)}^{- 1}{{diag}\left( Z_{S} \right)}}}{{\overset{\sim}{r}}_{x} = {r_{x} + {{{diag}(X)}^{- 1}r_{Z_{x}}}}}{{\overset{\sim}{r}}_{s} = {r_{s} + {{{diag}(S)}^{- 1}r_{Z_{s}}}}}{{\Delta\; X} = \begin{bmatrix} \left\{ {\Delta\; x_{i}} \right\} \\ \left\{ {\Delta\; x_{ij}} \right\} \\ \left\{ {\Delta\; u_{i}} \right\} \\ \left\{ {\Delta\; u_{ij}} \right\} \end{bmatrix}}{{\Delta\; y_{1}} = \left\{ {\Delta\; y_{1}^{x_{i}}} \right\}}{{\Delta\; y_{2}} = \left\{ {\Delta\; y_{2}^{x_{ij}}} \right\}}} & (29) \end{matrix}$ one can obtain, according to an embodiment of the invention:

$\begin{matrix} {{{\begin{bmatrix} H_{41} & 0 \\ H_{51} & 0 \\ {- I} & {- I} \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} = \begin{bmatrix} {- r_{y\; 1}} \\ {- r_{y\; 2}} \\ {- r_{Y}} \end{bmatrix}},} & (30) \\ {and} & \; \\ {{{\begin{bmatrix} {\overset{\sim}{H}}_{11} & 0 \\ 0 & {\overset{\sim}{H}}_{22} \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} + {\begin{bmatrix} {- I} & H_{41}^{T} & H_{51}^{T} \\ {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix}}} = {\begin{bmatrix} {- {\overset{\sim}{r}}_{x}} \\ {- {\overset{\sim}{r}}_{s}} \end{bmatrix}.}} & (31) \end{matrix}$ More specifically, EQS. (30) and (31) can be expanded using the definitions of EQS. (24.1)-(24.4) and the model definitions of EQS. (14.1)-(14.5) and (15.1)-(15.5). Thus, EQ. (30) has the following expanded form for each bus i and branch(i, j):

$\begin{matrix} {{{{\Delta\; x_{i}} - {\sum\limits_{j}^{\;}\;{R_{ij}\Delta\; x_{ij}}}} = r_{y_{1}}^{x_{i}}}{{{\Delta\; x_{ij}} - {R_{ij}\left( {{A_{ij}\Delta\; u_{i}} + {B_{ij}\Delta\; u_{j}} + {C_{ij}\Delta\; u_{ij}} + {D_{ij}\Delta\; u_{ji}}} \right)}} = r_{y_{2}}^{x_{ij}}}{{{\Delta\; X} + {\Delta\; S}} = r_{Y}}} & (32) \end{matrix}$ Similarly, EQ. (31) has the following expanded form:

$\begin{matrix} {\mspace{76mu}{{{{{\overset{\sim}{V}}_{i}^{X}\Delta\; x_{i}} - {\Delta\; Y^{x_{i}}} - {\Delta\; y_{1}^{x_{i}}}} = {- {\overset{\sim}{r}}_{x_{i}}}}\mspace{76mu}{{{{\overset{\sim}{W}}_{ij}^{X}\Delta\; x_{ij}} - {\Delta\; Y^{x_{ij}}} + {R_{ij}\Delta\; y_{1}^{x_{i}}} - {\Delta\; y_{2}^{x_{ij}}}} = {- {\overset{\sim}{r}}_{x_{ij}}}}{{{{\overset{\sim}{V}}_{i}^{U}\Delta\; u_{i}} - {\Delta\; Y^{u_{i}}} + {\sum\limits_{j}^{\;}\;{R_{ij}A_{ij}^{T}\Delta\; y_{2}^{x_{ij}}}} + {\sum\limits_{j}^{\;}\;{R_{ji}B_{ji}^{T}\Delta\; y_{2}^{x_{ji}}}}} = {- {\overset{\sim}{r}}_{u_{i}}}}\mspace{76mu}{{{{\overset{\sim}{W}}_{ij}^{U}\Delta\; u_{ij}} - {\Delta\; Y^{u_{ij}}} + {R_{ij}C_{ij}^{T}\Delta\; y_{2}^{x_{ij}}} + {R_{ji}D_{ji}^{T}\Delta\; y_{2}^{x_{ji}}}} = {- {\overset{\sim}{r}}_{u_{ij}}}}\mspace{76mu}{{{{\overset{\sim}{H}}_{22}\Delta\; S} - {\Delta\; Y}} = {- {\overset{\sim}{r}}_{S}}}}} & (33) \end{matrix}$

According to an embodiment of the invention, by first solving for control variables ΔU=[{Δu_(i)}, {Δu_(ij)}]^(T) in EQS. (33), the other variables in EQS. (32) and (33) can be solved using the following sequence of substitutions: {Δu _(i)},{Δu_(ij) }→{Δx _(ij) }→{Δx _(i) }→ΔS→ΔY→Δy ₁ →Δy ₂.  (34)

FIG. 2 is a flowchart of a method according to an embodiment of the invention for estimating the optimal power flow of a smart electric power grid. Referring now to the figure, a method according to an embodiment of the invention begins at step 21 by providing a cost function that models a smart electric power grid, such as the cost functions of EQS. (7) and (13), above, and EQ. (64), below. Constraints for the bus and branch control vectors, and the bus and branch dependent variables vectors, are provided at step 22, based on physical and operational requirements and actual measurements. At step 23, a set of linear equations is derived that minimize the cost function subject to the constraints. The dimension of the linear equations can be reduced at step 24 by solving for a subset of the linear equations. The reduced dimension linear equations are re-organized into primal and dual parts at step 25, and the re-organized reduced dimension linear equations are de-composed at step 26 into two systems of block matrix equations that can solved. A solution of the two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through the smart electric power grid.

Coding Details

According to an embodiment of the invention, a preconditioned conjugate gradient method can be used to solve the above equations, which uses a Jacobi preconditioning matrix composed of the diagonal elements of a matrix M of the linear system Mx=b. For every Δu_(i), the diagonal parameter can be computed as:

$M^{u_{i}} = {{\overset{\sim}{V}}_{i}^{U} + I + {\sum\limits_{j}^{\;}\;{R_{ij} \cdot A_{ij}^{T} \cdot \left\lbrack {{\left( {{\overset{\sim}{W}}_{ij}^{X} + {\overset{\sim}{H}}_{22}^{x_{ij}}} \right)A_{ij}} + {\left( {{\overset{\sim}{V}}_{i}^{X} + {\overset{\sim}{H}}_{22}^{x_{i}}} \right)\left( {\sum{R_{ij}A_{ij}}} \right)}} \right\rbrack}} + {\sum{{R_{ji} \cdot {B_{ji}^{T}\left( {{\overset{\sim}{W}}_{ji}^{X} + {\overset{\sim}{H}}_{22}^{x_{ji}} + {\overset{\sim}{V}}_{j}^{X} + {\overset{\sim}{H}}_{22}^{x_{j}}} \right)}}B_{ji}}}}$ Similarly, for every Δu_(ij), the diagonal parameter can be computed as:

$M^{u_{ij}} = {{\overset{\sim}{W}}_{ij}^{U} + {\overset{\sim}{H}}_{22}^{u_{ij}} + {{R_{ij} \cdot {C_{ij}^{T}\left( {{\overset{\sim}{W}}_{ij}^{X} + {\overset{\sim}{H}}_{22}^{x_{ij}} + {\overset{\sim}{V}}_{i}^{X} + {\overset{\sim}{H}}_{22}^{x_{i}}} \right)}}C_{ij}} + {{R_{ji} \cdot {D_{ji}^{T}\left( {{\overset{\sim}{W}}_{ji}^{X} + {\overset{\sim}{H}}_{22}^{x_{ji}} + {\overset{\sim}{V}}_{j}^{X} + {\overset{\sim}{H}}_{22}^{x_{j}}} \right)}}D_{ji}^{T}}}$

According to an embodiment of the invention, the quadratic objective function ƒ(X) is written as ƒ(X)=½X^(T)QX+p^(T)X, which has the following associated dual program and duality gap:

max   P^(T)Q⁻¹P^(T) + y₁^(T)b_(x) + y₂^(T)b_(v) + Y^(T)X^(UB) $s.t.\mspace{14mu}\left\{ \begin{matrix} {{Z_{S} + Y} \leq 0} \\ {{Q^{- 1}P} \geq 0} \\ {Z_{S} \geq 0} \\ {Z_{X} \geq 0} \end{matrix} \right.$ where P=p−Z _(x) −Y−y ₁ ^(T) ∇o(X)−y ₂ ^(T) ∇v(X). According to an embodiment of the invention, the current duality gap is: ƒ(X)−g(Y,u ₁ ,y ₂ ,Z _(S) ,Z _(X))=Y·(−r _(Y))+y ₁(−r _(y) ₁ )+y ₂·(−r _(y) ₂ )+½r _(X) ^(T) Q ⁻¹ r _(X) ^(T)+(Z _(X) ·X+Z _(S) ·S)−S·(Z _(S) +Y). In this expression of the duality gap, the first four terms come from the primal infeasibility for the current primal variable (X, S), and the last two terms come from the complementary slackness. When an optimal primal and dual solution is reached, the duality gap is 2nμ, where n is the dimension. Interior Point Method Solver Math for a Discrete AC OPF Model

First, a model according to an embodiment of the invention is standardized to the following format:

$\begin{matrix} {\mspace{79mu}{{X + S_{X} - X^{UB}} = {0\mspace{14mu}\left\{ \begin{matrix} {{x_{i} + s_{x_{i}}} = x_{i}^{UB}} \\ {{x_{ij} + s_{x_{ij}}} = x_{ij}^{UB}} \end{matrix} \right.}}} & (35.1) \\ {\mspace{79mu}{{{MU} + {NZ} + S_{U}} = {b\mspace{14mu}\left\{ \begin{matrix} {{{\sum\limits_{Ki}u_{i}^{k}} - s_{ui}^{1,k}} = u_{i}^{LB}} \\ {{{\sum\limits_{Ki}u_{i}^{k}} + s_{ui}^{2,k}} = u_{i}^{UB}} \\ {{{\sum\limits_{Ki}u_{ij}^{k}} - s_{uij}^{1,k}} = u_{ij}^{LB}} \\ {{{\sum\limits_{Ki}u_{ij}^{k}} + s_{uij}^{1,k}} = u_{ij}^{UB}} \\ {{{\sum\limits_{Ki}z_{i}^{k}} - s_{zi}^{1,k}} = z_{i}^{LB}} \\ {{{\sum\limits_{Ki}z_{i}^{k}} + s_{zi}^{2,k}} = z_{i}^{UB}} \\ {{{\sum\limits_{Ki}z_{ij}^{k}} - s_{zij}^{1,k}} = z_{ij}^{LB}} \\ {{{\sum\limits_{Ki}z_{ij}^{k}} + s_{zij}^{2,k}} = z_{ij}^{UB}} \\ {{{- u_{i}^{k}} + {z_{i}^{k} \cdot u_{ik}^{LB}} + s_{uzi}^{1,k}} = 0} \\ {{u_{i}^{k} - {z_{i}^{k} \cdot u_{ik}^{UB}} + s_{uzi}^{2,k}} = 0} \\ {{{- u_{ij}^{k}} + {z_{ij}^{k} \cdot u_{ijk}^{LB}} + s_{uzij}^{1,k}} = 0} \\ {{u_{i}^{k} - {z_{ij}^{k} \cdot u_{ijk}^{UB}} + s_{uzij}^{2,k}} = 0} \end{matrix} \right.}}} & (35.2) \\ {\mspace{79mu}{{{{\alpha(X)}\text{:}\mspace{14mu} x_{i}} - {\sum\limits_{j}^{\;}{R_{ij}x_{ij}}} - b_{v}} = 0}} & (35.3) \\ {{{{\beta\left( {X,U} \right)}\text{:}\mspace{14mu} x_{ij}} - {R_{ij} \cdot \left( {{\sum\limits_{K^{i}}^{\;}{A_{ij}^{k} \cdot u_{i}^{k}}} + {\sum\limits_{K^{i}}^{\;}{B_{ij}^{k} \cdot u_{j}^{k}}} + {\sum\limits_{K^{ij}}^{\;}{C_{ij}^{k} \cdot u_{ij}^{k}}} + {\sum\limits_{K^{ji}}^{\;}{D_{ij}^{k} \cdot u_{ji}^{k}}}} \right)}} = 0} & (35.4) \\ {\mspace{79mu}{{{{\gamma\left( {X,U,Z} \right)}\text{:}\mspace{14mu}{FX}} + {GU} + {HZ} + J} = 0}} & (35.5) \\ {\mspace{79mu}\left( {{Cutting}\mspace{14mu}{planes}} \right)} & \; \end{matrix}$ where the variables X, U, Z, M, N, and b are defined as follows:

$\begin{matrix} \begin{matrix} {X = \begin{bmatrix} \left\{ x_{i} \right\} \\ \left\{ x_{ij} \right\} \end{bmatrix}} & \; \\ {U = \begin{bmatrix} \left\{ u_{i}^{k} \right\} \\ \left\{ x_{ij}^{k} \right\} \end{bmatrix}} & \; \\ {Z = \begin{bmatrix} \left\{ z_{i}^{k} \right\} \\ \left\{ z_{ij}^{k} \right\} \end{bmatrix}} & \; \end{matrix} & (36.1) \\ {M = \begin{bmatrix} {- E_{1}} & 0 \\ E_{1} & 0 \\ 0 & {- E_{1}} \\ 0 & E_{1} \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ {- I} & 0 \\ I & 0 \\ 0 & {- I} \\ 0 & I \end{bmatrix}} & (36.2) \\ {N = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ {- E_{1}} & 0 \\ E_{1} & 0 \\ 0 & {- E_{1}} \\ 0 & E_{1} \\ I_{1} & 0 \\ {- I_{2}} & 0 \\ 0 & I_{3} \\ 0 & I_{4} \end{bmatrix}} & \; \\ {b = \begin{bmatrix} \left\{ {- u_{i}^{LB}} \right\} \\ \left\{ u_{i}^{UB} \right\} \\ \left\{ {- u_{ij}^{LB}} \right\} \\ \left\{ u_{ij}^{UB} \right\} \\ \left\{ {- z_{i}^{LB}} \right\} \\ \left\{ z_{i}^{UB} \right\} \\ \left\{ {- z_{ij}^{LB}} \right\} \\ \left\{ z_{ij}^{UB} \right\} \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}} & \; \end{matrix}$ and where

$\begin{matrix} {E_{1} = \begin{bmatrix} {\mathbb{e}}^{T} & 0 & \ldots & 0 \\ 0 & {\mathbb{e}}^{T} & \ldots & 0 \\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & \ldots & {\mathbb{e}}^{T} \end{bmatrix}} & (36.3) \\ {I_{1} = {{diag}\left\{ u_{ik}^{LB} \right\}}} & (36.4) \\ {I_{2} = {{diag}\left\{ u_{ik}^{UB} \right\}}} & \; \\ {I_{3} = {{diag}\left\{ u_{ijk}^{LB} \right\}}} & \; \\ {I_{4} = {{diag}\left\{ u_{ijk}^{uB} \right\}}} & \; \end{matrix}$

Following the above notation, a Lagrangian function according to an embodiment of the invention can be formulated as follows: L=ƒ(X,U)+μ^(T)(ln X+ln U+ln Z+ln S _(X)+ln S _(U))+λ₁ ^(T)(X+S _(X) −X ^(UB))+λ₂ ^(T)(MU+NZ+S _(U) −b)+y ₁ ^(T)α(X)+y ₂ ^(T)β(X,U)+y ₃ ^(T)γ(X,U,Z)  (37) The KKT first order optimality conditions of the resulting program are:

$\begin{matrix} {\begin{bmatrix} {\nabla L_{X}} \\ {\nabla L_{U}} \\ {\nabla L_{Z}} \\ {\nabla L_{S_{X}}} \\ {\nabla L_{S_{U}}} \\ {\nabla L_{\lambda_{1}}} \\ {\nabla L_{\lambda_{2}}} \\ {\nabla L_{y_{1}}} \\ {\nabla L_{y_{2}}} \\ {\nabla L_{y_{3}}} \end{bmatrix} = {\quad{\quad{\begin{bmatrix} {{f_{X}\left( {X,U} \right)} + {\mu/X} + \lambda_{1}^{T} + {y_{1}^{T}{\alpha_{X}(X)}} + {y_{2}^{T}{\beta_{X}\left( {X,U} \right)}} +} \\ {y_{3}^{T}{\gamma_{X}\left( {X,U,Z} \right)}} \\ {{f_{U}\left( {X,U} \right)} + {\mu/U} + {y_{2}^{T}{\beta_{U}\left( {X,U} \right)}} + {y_{3}^{T}{\gamma_{U}\left( {X,U,Z} \right)}}} \\ {{\mu/Z} + {\lambda_{2}^{T}N} + {y_{3}^{T}{\gamma_{Z}\left( {X,U,Z} \right)}}} \\ {{\mu/S_{X}} + \lambda_{1}^{T}} \\ {{\mu/S_{U}} + \lambda_{2}^{T}} \\ {X + S_{X} - X^{UB}} \\ {{MU} + {NZ} + S_{U} - b} \\ {\alpha(X)} \\ {\beta\left( {X,U} \right)} \\ {\gamma\left( {X,U,Z} \right)} \end{bmatrix}{\quad{= 0}}}}}} & (38) \end{matrix}$ By introducing

$\begin{matrix} {T = {\begin{bmatrix} T_{X} \\ T_{U} \\ T_{Z} \\ T_{S_{X}} \\ T_{S_{U}} \end{bmatrix} = \begin{bmatrix} {\mu/X} \\ {\mu/U} \\ {\mu/Z} \\ {\mu/S_{X}} \\ {\mu/S_{U}} \end{bmatrix}}} & (39) \end{matrix}$ The optimality conditions can be rewritten as

$\begin{matrix} {\begin{bmatrix} \begin{matrix} {{f_{X}\left( {X,U} \right)} + T_{X} + \lambda_{1}^{T} + {y_{1}^{T}\alpha_{X}(X)} +} \\ {{y_{2}^{T}{\beta_{X}\left( {X,U} \right)}} + {y_{3}^{T}{\gamma_{X}\left( {X,U,Z} \right)}}} \end{matrix} \\ {{f_{U}\left( {X,U} \right)} + T_{U} + {y_{2}^{T}{\beta_{U}\left( {X,U} \right)}} + {y_{3}^{T}{\gamma_{U}\left( {X,U,Z} \right)}}} \\ {T_{Z} + {\lambda_{2}^{T}N} + {y_{3}^{T}{\gamma_{Z}\left( {X,U,Z} \right)}}} \\ {T_{S_{X}} + \lambda_{1}^{T}} \\ {T_{S_{U}} + \lambda_{2}^{T}} \\ {X + S_{X} - X^{UB}} \\ {{MU} + {NZ} + S_{U} - b} \\ {\alpha(X)} \\ {\beta\left( {X,U} \right)} \\ {\gamma\left( {X,U,Z} \right)} \\ {{T \cdot \begin{bmatrix} X \\ U \\ Z \\ S_{X} \\ S_{U} \end{bmatrix}} - {\mu\; e}} \end{bmatrix} = 0.} & (40) \end{matrix}$ Newton's method can be used to solve the above non-linear equation, to yield the following linear equation:

$\begin{matrix} {{H \cdot \begin{bmatrix} {\Delta\; X} \\ {\Delta\; U} \\ {\Delta\; Z} \\ {\Delta\; S_{X}} \\ {\Delta\; S_{U}} \\ {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; y_{3}} \\ {\Delta\; T} \end{bmatrix}} = {- \begin{bmatrix} \begin{matrix} {{f_{X}\left( {X,U} \right)} + T_{X} + \lambda_{1}^{T} + {y_{1}^{T}\alpha_{X}(X)} +} \\ {{y_{2}^{T}{\beta_{X}\left( {X,U} \right)}} + {y_{3}^{T}{\gamma_{X}\left( {X,U,Z} \right)}}} \end{matrix} \\ {{f_{U}\left( {X,U} \right)} + T_{U} + {y_{2}^{T}{\beta_{U}\left( {X,U} \right)}} + {y_{3}^{T}{\gamma_{U}\left( {X,U,Z} \right)}}} \\ {T_{Z} + {\lambda_{2}^{T}N} + {y_{3}^{T}{\gamma_{Z}\left( {X,U,Z} \right)}}} \\ {T_{S_{X}} + \lambda_{1}^{T}} \\ {T_{S_{U}} + \lambda_{2}^{T}} \\ {X + S_{X} - X^{UB}} \\ {{MU} + {NZ} + S_{U} - b} \\ {\alpha(X)} \\ {\beta\left( {X,U} \right)} \\ {\gamma\left( {X,U,Z} \right)} \\ {{T \cdot \begin{bmatrix} X \\ U \\ Z \\ S_{X} \\ S_{U} \end{bmatrix}} - {\mu\; e}} \end{bmatrix}}} & (41) \end{matrix}$ where H is the following Jacobian matrix:

$\begin{matrix} {H = \begin{bmatrix} H_{1,1} & 0 & 0 & 0 & 0 & I & 0 & H_{1,8} & H_{1,9} & H_{1,10} & E_{1} \\ 0 & H_{2,2} & 0 & 0 & 0 & 0 & H_{2,7} & 0 & H_{2,9} & H_{2,10} & E_{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & H_{3,7} & 0 & 0 & H_{3,10} & E_{3} \\ 0 & 0 & 0 & 0 & 0 & I & 0 & 0 & 0 & 0 & E_{4} \\ 0 & 0 & 0 & 0 & 0 & 0 & I & 0 & 0 & 0 & E_{5} \\ I & 0 & 0 & I & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & H_{2,7} & H_{3,7} & 0 & I & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,8} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,9} & H_{2,9} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,10} & H_{2,10} & H_{3,10} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ T_{1} & T_{2} & T_{3} & T_{4} & T_{5} & 0 & 0 & 0 & 0 & 0 & T_{6} \end{bmatrix}} & (42) \end{matrix}$ where E_(i) is the i^(th) block identity matrix where

$\begin{matrix} {\mspace{79mu}{{I = {\sum\limits_{i = 1}^{5}E_{i}}};}} & (43.1) \\ {{T_{i} = {E_{i}\begin{bmatrix} {{diag}\left\{ T_{X} \right\}} & \; & \; & \; & \; \\ \; & {{diag}\left\{ T_{U} \right\}} & \; & \; & \; \\ \; & \; & {{diag}\left\{ T_{Z} \right\}} & \; & \; \\ \; & \; & \; & {{diag}\left\{ T_{S_{X}} \right\}} & \; \\ \; & \; & \; & \; & {{diag}\left\{ T_{S_{U}} \right\}} \end{bmatrix}}}\mspace{20mu}{{{{for}\mspace{14mu} i} = 1},2,3,4,{5;}}} & (43.2) \\ {\mspace{79mu}{{T_{6} = \begin{bmatrix} {{diag}\left\{ X \right\}} & \; & \; & \; & \; \\ \; & {{diag}\left\{ U \right\}} & \; & \; & \; \\ \; & \; & {{diag}\left\{ Z \right\}} & \; & \; \\ \; & \; & \; & {{diag}\left\{ S_{X} \right\}} & \; \\ \; & \; & \; & \; & {{diag}\left\{ S_{U} \right\}} \end{bmatrix}};}} & (43.3) \\ {\mspace{79mu}{{H_{1,1} = \left\lbrack V^{x} \middle| W^{x} \right\rbrack};}} & (43.4) \\ {\mspace{79mu}{{H_{1,8} = \left\lbrack I \middle| {- R} \right\rbrack^{T}};}} & (43.5) \\ {{R = \left( \begin{matrix} \left\lbrack {R_{11}I\mspace{14mu}\ldots\mspace{14mu} R_{1N_{1}}I} \right\rbrack & 0 & \ldots & 0 \\ 0 & \left\lbrack {R_{21}I\mspace{14mu}\ldots\mspace{14mu} R_{2N_{2}}I} \right\rbrack & \ldots & 0 \\ 0 & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & \left\lbrack {R_{N_{1}1}I\mspace{14mu}\ldots\mspace{14mu} R_{N_{1}N_{1}}I} \right\rbrack \end{matrix} \right)};} & (43.6) \\ {\mspace{79mu}{{H_{1.9} = \left\lbrack 0 \middle| I \right\rbrack^{T}};}} & (43.7) \\ {\mspace{79mu}{{H_{1.10} = F^{T}};}} & (43.8) \\ {\mspace{79mu}{{H_{2.2} = \left\lbrack V^{U} \middle| W^{U} \right\rbrack};}} & (43.9) \\ {\mspace{79mu}{{H_{2.7} = M^{T}};}} & (43.10) \\ {\mspace{79mu}{{H_{2.9} = \left\lbrack A \middle| B \right\rbrack^{T}},}} & (43.11) \end{matrix}$ where A_(ij), B_(ij) that correspond to the branch (i,j) are defined as A _(ij)=[0, . . . A _(ij) ¹ , . . . A _(ij) ^(K) ^(i) , 0, . . . B _(ij) ¹ , . . . B _(ij) ^(K) ^(j) ,0, . . . ],  (43.12) B _(ij)=[0, . . . C _(ij) ¹ , . . . C _(ij) ^(K) ^(ij) ,0, . . . D _(ij) ¹ , . . . D _(ij) ^(K) ^(ji) ,0, . . . ],  (43.13) H _(2.10) =N ^(T);  (43.14) H _(3.10) =H ^(T).  (43.15) By defining

$\begin{matrix} {\begin{bmatrix} r_{X} \\ r_{U} \\ r_{Z} \\ r_{S_{X}} \\ r_{S_{U}} \\ r_{\lambda_{1}} \\ r_{\lambda_{2}} \\ r_{y_{1}} \\ r_{y_{2}} \\ r_{y_{3}} \\ r_{T} \end{bmatrix} = {\quad\left\lbrack \begin{matrix} \begin{matrix} {{f_{X}\left( {X,U} \right)} + T_{X} + \lambda_{1}^{T} + {y_{1}^{T}\alpha_{X}(X)} +} \\ {{y_{2}^{T}{\beta_{X}\left( {X,U} \right)}} + {y_{3}^{T}{\gamma_{X}\left( {X,U,Z} \right)}}} \end{matrix} \\ {{f_{U}\left( {X,U} \right)} + T_{U} + {y_{2}^{T}{\beta_{U}\left( {X,U} \right)}} + {y_{3}^{T}{\gamma_{U}\left( {X,U,Z} \right)}}} \\ {T_{Z} + {\lambda_{2}^{T}N} + {y_{3}^{T}{\gamma_{Z}\left( {X,U,Z} \right)}}} \\ {T_{S_{X}} + \lambda_{1}^{T}} \\ {T_{S_{U}} + \lambda_{2}^{T}} \\ {X + S_{X} - X^{UB}} \\ {{MU} + {NZ} + S_{U} - b} \\ {\alpha(X)} \\ {\beta\left( {X,U} \right)} \\ {\gamma\left( {X,U,Z} \right)} \\ {{T \cdot \begin{bmatrix} X \\ U \\ Z \\ S_{X} \\ S_{U} \end{bmatrix}} - {\mu\; e}} \end{matrix} \right\rbrack}} & (44) \end{matrix}$ one has

$\begin{matrix} {\begin{bmatrix} H_{1,1} & 0 & 0 & 0 & 0 & I & 0 & H_{1,8} & H_{1,9} & H_{1,10} & E_{1} \\ 0 & H_{2,2} & 0 & 0 & 0 & 0 & H_{2,7} & 0 & H_{2,9} & H_{2,10} & E_{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & H_{3,7} & 0 & 0 & H_{3,10} & E_{3} \\ 0 & 0 & 0 & 0 & 0 & I & 0 & 0 & 0 & 0 & E_{4} \\ 0 & 0 & 0 & 0 & 0 & 0 & I & 0 & 0 & 0 & E_{5} \\ I & 0 & 0 & I & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & H_{2,7} & H_{3,7} & 0 & I & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,8} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,9} & H_{2,9} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,10} & H_{2,10} & H_{3,10} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ T_{1} & T_{2} & T_{3} & T_{4} & T_{5} & 0 & 0 & 0 & 0 & 0 & T_{6} \end{bmatrix}{\quad{\begin{bmatrix} {\Delta\; X} \\ {\Delta\; U} \\ {\Delta\; Z} \\ {\Delta\; S_{X}} \\ {\Delta\; S_{U}} \\ {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; y_{3}} \\ {\Delta\; T} \end{bmatrix} = {- {\begin{bmatrix} r_{X} \\ r_{U} \\ r_{Z} \\ r_{S_{X}} \\ r_{S_{U}} \\ r_{\lambda_{1}} \\ r_{\lambda_{2}} \\ r_{y_{1}} \\ r_{y_{2}} \\ r_{y_{3}} \\ r_{T} \end{bmatrix}.}}}}} & (45) \end{matrix}$

Similar to EQS. (26), the dimension of the system of EQS. (45) can be very large. By decomposing r_(T) as follows,

$\begin{matrix} {{r_{T} = \begin{bmatrix} r_{Tx} \\ r_{T_{U}} \\ r_{Tz} \\ r_{{Ts}_{X}} \\ r_{{Ts}_{U}} \end{bmatrix}},} & (46) \end{matrix}$ the dimension of EQS. (45) can be reduced as follows. From the last block constraints there is:

$\begin{matrix} {\begin{bmatrix} {\Delta\; T_{X}} \\ {\Delta T}_{U} \\ {\Delta\; T_{Z}} \\ {\Delta\; T_{S_{X}}} \\ {\Delta\; T_{S_{U}}} \end{bmatrix} = {\begin{bmatrix} {{{diag}^{- 1}(X)}\left( {r_{T_{X}} - {{{diag}\left( T_{X} \right)}\Delta\; X}} \right)} \\ {{{diag}^{- 1}(U)}\left( {r_{T_{U}} - {{{diag}\left( T_{U} \right)}\Delta\; U}} \right)} \\ {{{diag}^{- 1}(Z)}\left( {r_{T_{Z}} - {{{diag}\left( T_{Z} \right)}\Delta\; Z}} \right)} \\ {{{diag}^{- 1}\left( S_{X} \right)}\left( {r_{T_{S_{X}}} - {{{diag}\left( T_{S_{X}} \right)}\Delta\; S_{X}}} \right)} \\ {{{diag}^{- 1}\left( S_{U} \right)}\left( {r_{T_{S_{U}}} - {{{diag}\left( T_{S_{U}} \right)}\Delta\; S_{U}}} \right)} \end{bmatrix}.}} & (47) \end{matrix}$ This definition for the ΔTs can be substituted into the other blocks of the matrix on the left hand side of EQ. (45), which can be reorganized into primal and dual parts:

$\begin{matrix} {\begin{bmatrix} I & 0 & 0 & I & 0 & ❘ & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & H_{2,7} & H_{3,7} & 0 & I & ❘ & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,8} & 0 & 0 & 0 & 0 & ❘ & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,9} & H_{2,9} & 0 & 0 & 0 & ❘ & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{1,10} & H_{2,10} & H_{3,10} & 0 & 0 & ❘ & 0 & 0 & 0 & 0 & 0 & 0 \\  - & - & - & - & - & - & - & - & - & - & - & \; \\ {\overset{\sim}{H}}_{1,1} & 0 & 0 & 0 & 0 & ❘ & I & 0 & H_{1,8} & H_{1,9} & H_{1,10} & E_{1} \\ 0 & {\overset{\sim}{H}}_{2,2} & 0 & 0 & 0 & ❘ & 0 & H_{2,7} & 0 & H_{2,9} & H_{2,10} & E_{2} \\ 0 & 0 & {\overset{\sim}{H}}_{3,3} & 0 & 0 & ❘ & 0 & H_{3,7} & 0 & 0 & H_{3,10} & E_{3} \\ 0 & 0 & 0 & {\overset{\sim}{H}}_{4,4} & 0 & ❘ & I & 0 & 0 & 0 & 0 & E_{4} \\ 0 & 0 & 0 & 0 & {\overset{\sim}{H}}_{5,5} & ❘ & 0 & I & 0 & 0 & 0 & E_{5} \end{bmatrix}{\quad{\begin{bmatrix} {\Delta\; X} \\ {\Delta\; U} \\ {\Delta\; Z} \\ {\Delta\; S_{X}} \\ {\Delta\; S_{U}} \\ {- ❘} \\ {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; y_{3}} \end{bmatrix} = {{- \begin{bmatrix} r_{\lambda_{1}} \\ r_{\lambda_{2}} \\ r_{y_{1}} \\ r_{y_{2}} \\ r_{y_{3}} \\  - \\ {\overset{\sim}{r}}_{X} \\ {\overset{\sim}{r}}_{U} \\ {\overset{\sim}{r}}_{Z} \\ {\overset{\sim}{r}}_{S_{X}} \\ {\overset{\sim}{r}}_{S_{U}} \end{bmatrix}}\mspace{20mu}{where}}}}} & (48) \\ {\mspace{79mu}{{\begin{bmatrix} {\overset{\sim}{H}}_{1,1} \\ {\overset{\sim}{H}}_{2,2} \\ {\overset{\sim}{H}}_{3,3} \\ {\overset{\sim}{H}}_{4,4} \\ {\overset{\sim}{H}}_{5,5} \end{bmatrix} = \begin{bmatrix} {H_{1,1} - {{{diag}^{- 1}(X)}{{diag}\left( T_{X} \right)}}} \\ {H_{2,2} - {{{diag}^{- 1}(U)}{{diag}\left( T_{U} \right)}}} \\ {{- {{diag}^{- 1}(Z)}}{{diag}\left( T_{Z} \right)}} \\ {{- {{diag}^{- 1}\left( S_{X} \right)}}{{diag}\left( T_{S_{X}} \right)}} \\ {{- {{diag}^{- 1}\left( S_{U} \right)}}{{diag}\left( T_{S_{U}} \right)}} \end{bmatrix}},\mspace{20mu}{and}}} & (49) \\ {\mspace{79mu}{\begin{bmatrix} {\overset{\sim}{r}}_{X} \\ {\overset{\sim}{r}}_{U} \\ {\overset{\sim}{r}}_{Z} \\ {\overset{\sim}{r}}_{S_{X}} \\ {\overset{\sim}{r}}_{S_{U}} \end{bmatrix} = {\begin{bmatrix} {r_{X} + {{{diag}^{- 1}(X)}r_{T_{X}}}} \\ {r_{U} + {{{diag}^{- 1}(U)}r_{T_{U}}}} \\ {r_{Z} + {{{diag}^{- 1}(Z)}r_{T_{Z}}}} \\ {r_{S_{X}} + {{{diag}^{- 1}\left( S_{X} \right)}r_{S_{T_{X}}}}} \\ {r_{S_{U}} + {{{diag}^{- 1}\left( S_{U} \right)}r_{S_{T_{U}}}}} \end{bmatrix}.}}} & (50) \end{matrix}$ According to an embodiment of the invention, EQ. (45) can be reduced to

$\begin{matrix} {\mspace{79mu}{{{\begin{bmatrix} I & 0 & 0 & I & 0 \\ 0 & H_{2,7} & H_{3,7} & 0 & I \\ H_{1,8} & 0 & 0 & 0 & 0 \\ H_{1,9} & H_{2,9} & 0 & 0 & 0 \\ H_{1,10} & H_{2,10} & H_{3,10} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; U} \\ {\Delta\; Z} \\ {\Delta\; S_{X}} \\ {\Delta\; S_{U}} \end{bmatrix}} = {- \begin{bmatrix} r_{\lambda_{1}} \\ r_{\lambda_{2}} \\ r_{y_{1}} \\ r_{y_{2}} \\ r_{y_{3}} \end{bmatrix}}}\mspace{20mu}{and}}} & (51) \\ {\begin{bmatrix} {\overset{\sim}{H}}_{1,1} & 0 & 0 & 0 & 0 \\ 0 & {\overset{\sim}{H}}_{2,2} & 0 & 0 & 0 \\ 0 & 0 & {\overset{\sim}{H}}_{3,3} & 0 & 0 \\ 0 & 0 & 0 & {\overset{\sim}{H}}_{4,4} & 0 \\ 0 & 0 & 0 & 0 & {\overset{\sim}{H}}_{5,5} \end{bmatrix}{\quad{{\begin{bmatrix} {\Delta\; X} \\ {\Delta\; U} \\ {\Delta\; Z} \\ {\Delta\; S_{X}} \\ {\Delta\; S_{U}} \end{bmatrix} + {\begin{bmatrix} I & 0 & H_{1,8} & H_{1,9} & H_{1,10} & E_{1} \\ 0 & H_{2,7} & 0 & H_{2,9} & H_{2,10} & E_{2} \\ 0 & H_{3,7} & 0 & 0 & H_{3,10} & E_{3} \\ I & 0 & 0 & 0 & 0 & E_{4} \\ 0 & I & 0 & 0 & 0 & E_{5} \end{bmatrix}\begin{bmatrix} {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; y_{3}} \end{bmatrix}}} = \begin{bmatrix} {\overset{\sim}{r}}_{X} \\ {\overset{\sim}{r}}_{U} \\ {\overset{\sim}{r}}_{Z} \\ {\overset{\sim}{r}}_{S_{X}} \\ {\overset{\sim}{r}}_{S_{U}} \end{bmatrix}}}} & (52) \end{matrix}$ EQS. (51) and (52) can be solved similarly to EQS. (30) and (31) above, by expanding the block components in the matrices of EQS. (51) and (52) using their definitions above, and carrying out a series of substitutions similar to that of EQ. (34). Heuristic Block-Wise Newton Method

It can be seen that a system according to an embodiment of the invention is highly structural in the way that all constraints can be grouped by buses and branches. Without loss of generality, a minimal AC OPF model according to an embodiment of the invention will be considered as an example, and the same heuristic can be applied to a more general discrete AC OPF model according to an embodiment of the invention.

For a minimal AC OPF model according to an embodiment of the invention, the following non-linear system will be solved:

$\begin{matrix} {\overset{\sim}{L} = {\begin{bmatrix} {{\nabla{f(X)}} - {\mu/X} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{{- \mu}/S} - Y} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \end{bmatrix} = 0.}} & (53) \end{matrix}$ According to an embodiment of the invention, the vector {tilde over (L)} can be re-organized and rewritten as follows:

$\begin{matrix} {{\overset{\sim}{L} = \begin{bmatrix} \left\{ L^{i} \right\} \\ \left\{ L^{ij} \right\} \end{bmatrix}},} & (54) \end{matrix}$ where each L^(i) denotes the constraints of bus i, thus is a function of x_(i) and {x_(ik)}. Similarly each L^(ij) denotes the constraints of branch (i, j), which is a function of x_(i), x_(j), u_(i), u_(j), u_(ij), and u_(ji). Therefore, according to an embodiment of the invention, the system can be decomposed by different buses and branches. More specifically, instead of directly applying Newton method to {tilde over (L)}, it can be used for every bus L^(i) and branch L^(ij):

$\begin{matrix} {{{{H^{i} \cdot \begin{bmatrix} {\Delta\; X^{i}} \\ {\Delta\; S^{i}} \\ {\Delta\; Y^{i}} \\ {\Delta\; y_{1}^{i}} \\ {\Delta\; y_{2}^{i}} \end{bmatrix}} = {- L^{i}}};{{H^{ij} \cdot \begin{bmatrix} {\Delta\; X^{ij}} \\ {\Delta\; S^{ij}} \\ {\Delta\; Y^{ij}} \\ {\Delta\; y_{1}^{ij}} \\ {\Delta\; y_{2}^{ij}} \end{bmatrix}} = {- L^{ij}}}},} & (55) \end{matrix}$ where the Jacobian matrices H^(i) and H^(ij) can be obtained by the above techniques. Moreover, cutting plane methods can be used to solve these block-wise Newton's equations since each cutting plane is a function of either L^(i) or L^(ij).

One issue with respect to this block-wise Newton method according to an embodiment of the invention is its convergence because some variables will be updated more than once in each iteration. For instance, there will be different Δu_(i) when solving L^(ij) and L^(ik). A straightforward technique would be to add Au, whenever there is an update for it. However, no theoretic proof will guarantee the convergence in this case.

Interior Point Method Algorithm

Based on the above linear system, given a tolerance ε>0, barrier parameter μ=μ⁰, and iteration counter k=1, an interior point algorithm according to an embodiment of the invention can be designed as follows, with reference to the flowchart of FIG. 3.

Begin Initialize variables X^(k)=X⁰, S^(k)=S⁰, . . . (Step 31)

Repeat

-   -   Update μ^(k+1)←ε·μu^(k), εε(0.1, . . . , 0.9); (Step 32)     -   Evaluate (ΔX^(k), ΔS^(k), ΔY^(k), Δy_(i) ^(k), Δy₂ ^(k)) using         EQS. (30) and (31); (Step 33)     -   Calculate (ΔZ_(X) ^(k), ΔZ_(S) ^(k)) using EQS. (27); (Step 34)     -   Find a step length Δβ^(k) by a line search; (Step 35)     -   Update: (Step 36)         X ^(k+1) ←X ^(k)+Δβ^(k) ΔX ^(k),         S ^(k+1) ←S ^(k)+Δβ^(k) ΔS ^(k),         Y ^(k+1) ←Y ^(k)+Δβ^(k) ΔY ^(k),         y ₁ ^(k+1) ←y ₁ ^(k)+Δβ^(k) Δy ₁ ^(k),         y ₂ ^(k+1) ←y ₂ ^(k)+Δβ^(k) Δy ₂ ^(k),         Z _(X) ^(k+1) ←Z _(X) ^(k)+Δβ^(k) ΔZ _(X) ^(k),         Z _(S) ^(k+1) ←Z _(S) ^(k)+Δβ^(k) ΔZ _(S) ^(k),

Until Stop conditions are satisfied (Step 37).

The initial values of the variables depend upon the situation. To measure optimality, a dual gap according to an embodiment of the invention is defined as: gap=X ^(UB) ^(T) Y+b _(x) ^(T) y ₁ +b _(v) ^(T) y ₂−½X ^(T) VX−ƒ(x,u).  (56) A stop condition according to an embodiment of the invention can be one of the following: gap≦ε, or  (57) max{r _(X) ,r _(S)}≦ε.  (58)

The line search of Step 15 tries to find the suitable step size. Two methods which could be considered are the Trust-region method and the Filter line-search.

A trust region method has a region around the current search point, where the quadratic model q_(k)(β^(k))=ƒ(X_(k))+∇ƒ(X_(k))^(T)β^(k)+½β_(k) ^(T)B_(k)β^(k) for “local minimization” is “trusted” to be correct and steps are chosen to stay within this region. The size of the region is modified during the search, based on how well the model agrees with actual function evaluations. Once a step β^(k) is chosen, the function is evaluated at the new point, and the actual function value is checked against the value predicted by the quadratic model. What is actually computed is the ratio of actual to predicted reduction:

$\beta^{k} = {\frac{{f\left( X_{k} \right)} - {f\left( {X_{k} + \beta^{k}} \right)}}{{q_{k}(0)} - {q_{k}\left( \beta^{k} \right)}} = {\frac{{actual}\mspace{14mu}{reduction}\mspace{14mu}{of}\mspace{14mu} f}{{predicted}\mspace{14mu}{model}\mspace{14mu}{reduction}\mspace{14mu}{of}\mspace{14mu} f}.}}$ If β^(k) is close to 1, then the quadratic model is quite a good predictor and the region can be increased in size. On the other hand, if β^(k) is too small, the region is decreased in size. When k is below a threshold η, the step is rejected and recomputed.

A filter line-search interprets EQS. (19) and (14.1)-(14.5) as a bi-objective optimization program with twin goals of minimizing the objective function L_(μ)(X, S) and the constraint violation θ=∥c(X)∥. Following this paradigm, a trial point X^(k)(β^(k))=X^(k)+β^(k)ΔX^(k) may be acceptable, if it leads to sufficient progress toward either objective compared to the current iterate, i.e., if θ(X ^(k)(β^(k)))≦(1−γ_(θ))θ(X ^(k)), or L _(μ)(X ^(k)(β^(k)),S ^(k)(β^(k)))≦L _(μ)(X ^(k) ,S ^(k))−γ_(L)θ(X ^(k)) holds for fixed constant γ_(θ), γ_(L)ε(0,1). Extended OPF Model

An extended OPF model according to other embodiments of the invention assumes an AC-only grid with continuous-only dependent (non-control) variables x and control variables u, and that a particular bus can be connected to more than one generator or load, or both at the same time. In addition to that, no thermal capacity limits are enforced on the branch power flows. An extended OPF model according to other embodiments of the invention can be solved by the above disclosed methods for solving a discrete AC OPF model and a minimal AC OPF model according to other embodiments of the invention.

Bus dependent variables vector x_(i) typically includes net active x_(i) ^(P) and reactive x_(i) ^(Q), and it is computed for each bus i according to EQ. (1):

$\begin{matrix} {{x_{i} = {{\sum\limits_{k = 1}^{G_{i}}x_{i,k}^{g}} - {\sum\limits_{k = 1}^{L_{i}}x_{i,k}^{l}} + {\sum\limits_{k = 1}^{{SH}_{i}}x_{i,k}^{sh}}}},} & (60) \end{matrix}$ where G_(i), L_(i) and SH_(i) are the number of generators, loads and shunts placed on bus i and x_(i,k) ^(g), x_(i,k) ^(l) and x_(i,k) ^(sh) are the k^(th) active and reactive generation, load and shunt variables at bus i, respectively. Power system models may include loads and shunts as a function of the associated node's voltage magnitudes. However, in a model according to an embodiment of the invention, load and shunt devices are assumed to operate as constant power devices. EQ. (60) may be rewritten in a vector notation as x=Ã{tilde over (x)},  (61) where Ãε{−1,0,1}^(N×2(G+L+SH)) indicates which generators, loads or shunts are connected to the same bus, and {tilde over (x)} is a vector of generators, loads and shunts on the network.

Box Constraints:

Two box-constraint families may be defined that independently limit bus control vector u and dependent bus variable column vector {tilde over (x)}: {tilde over (x)} ^(LB) ≦{tilde over (x)}≦{tilde over (x)} ^(UB), u ^(LB) ≦u≦u ^(UB), Bus control vector u_(i) typically includes voltage magnitude V_(i) (log-scaled according to an embodiment of the invention), phase angle θ_(i); and shunt switch position Q_(i).

Similarly, box-constraints may be defined for branch variables that independently limit branch control vector v_(ij), which typically includes transformer tap selection t_(ij) and phase shift φ_(ij), and dependent branch variable vector y_(ij), which typically includes net active y_(ij) ^(P) and reactive y_(ij) ^(Q) power injections, for each branch ij: y ^(LB) ≦y≦y ^(UB), v ^(LB) ≦v≦v ^(UB),

Flow Balance Set:

At every bus, a net power injection balance should be maintained between a bus and the connected branches. A branch dependent variables vector y_(ij) may be defined, in which typically y_(ij)=y_(ij) ^(P), y_(ij) ^(Q)), and the balance for each bus i may be specified as follows:

$\begin{matrix} {{a_{i}\overset{\sim}{x}} = {\sum\limits_{j = 1}^{N}{R_{ij}y_{ij}}}} & (62) \end{matrix}$ where a_(i) is i^(th) row of matrix Ã. Note that matrix R is always highly sparse since the connection density in real-world power grids is typically 1 to 5 branches per bus.

Linearized Voltage Law:

According to an embodiment of the invention, a first order approximation of the voltage law constraints is used. A linearized version of the AC voltage law defining the relation between branch power flow vectors y_(ij) and control vectors u_(i), u_(j), v_(ij) which are linearized with respect to a given fixed point (u_(i) ⁰, u_(j) ⁰, v_(ij) ⁰) is written as y _(ij) =R _(ij)·(A _(ij) ⁰ ·u _(i) +B _(ij) ⁰ ·u _(i) +C _(ij) ⁰ ·v _(ij) +D _(ij) ⁰ ·v _(ji) +E _(ij) ⁰)  (63)

This constraint exhibits two types of sparsity: one inherited from R_(ij), and the other due to dependence on only two control vectors (u_(i), u_(j)) and (v_(ij), v_(ji)) instead of the whole-grid control vector U=└(u_(i))_(i=1) ^(N); (v_(ij))_(i,j=1) ^(N)┘. The matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i,j) for the given point (u_(i) ⁰, u_(j) ⁰, v_(ij) ⁰).

Objective Function:

An OPF model according to an embodiment of the invention assumes a generic convex quadratic cost function to be minimized. This function includes terms related to managing the vector u_(i) that minimizes control adjustment costs, bus dependent variables vector {tilde over (x)}_(i) that minimize generation costs, and branch dependent variables vector y_(ij) that minimize transmission losses, and a penalty coefficient for the auxiliary variable β: ƒ({tilde over (x)},u,y,v)=½[{tilde over (x)} ^(T) V ^(X) {tilde over (x)}+u ^(T) V ^(U) u+y ^(T) W ^(Y) y+v ^(T) W ^(V) v]+S ^(X) {tilde over (x)}+S ^(U) u+T ^(Y) y+T ^(V) v.  (64) Interior Point Method Solver Math for the Minimal OPF Model

A model according to an embodiment of the invention can be standardized to the following format: ½[x ₁ ^(T) V ^(X) {tilde over (x)} ₁ +u ₁ ^(T) V ^(U) u ₁ +y _(i) ^(T) W ^(Y) y ₁ +v ₁ ^(T) W ^(V) v ₁ ]+S ^(X) {tilde over (x)} ₁ +S ^(U) u ₁ +T ^(Y) y ₁ +T ^(V) v ₁.  (65) Slack variables are introduced for all bound constraints:

$\begin{matrix} {{{{\overset{\sim}{x}}_{i\; 1} + s_{x\; i}} = {{\overset{\sim}{x}}_{i\; 1}^{UB}{\forall i}}}{{u_{i\; 1} + s_{ui}} = {u_{i\; 1}^{UB}{\forall i}}}{{y_{{ij}\; 1} + s_{yij}} = {y_{{ij}\; 1}^{UB}{\forall{ij}}}}{{v_{{ij}\; 1} + s_{vij}} = {v_{{ij}\; 1}^{UB}{\forall{ij}}}}} & (66.1) \\ {{{{\overset{\sim}{\alpha}}_{i}{\overset{\sim}{x}}_{1}} - {\sum\limits_{j = 1}^{N}{R_{ij} \cdot y_{{ij}\; 1}}}} = b_{x}} & (66.2) \\ {{y_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot v_{{ij}\; 1}} + {D_{ij}^{0} \cdot v_{{ij}\; 1}}} \right)}} = b_{v}} & (66.3) \\ {{\overset{\sim}{x}}_{1},y_{1},u_{1},v_{1},s_{\overset{\sim}{x}},s_{u},s_{y},{s_{v} \geq 0}} & (66.4) \end{matrix}$ where S_({tilde over (x)}), S_(u), S_(y), S_(v) are slack variables, and

$\begin{matrix} {\mspace{20mu}{{{\overset{\sim}{x}}_{i\; 1} = {{\overset{\sim}{x}}_{i} - {\overset{\sim}{x}}_{i}^{LB}}},\mspace{20mu}{y_{{ij}\; 1} = {y_{ij} - y_{ij}^{LB}}},}} & (67.1) \\ {\mspace{20mu}{{{\overset{\sim}{x}}_{i\; 1}^{UB} = {{\overset{\sim}{x}}_{i}^{UB} - {\overset{\sim}{x}}_{i}^{LB}}},\mspace{20mu}{y_{i\; j\; 1}^{UB} = {y_{ij}^{UB} - y_{ij}^{LB}}},\mspace{20mu}{u_{i\; 1} = {u_{i} - u_{i}^{LB}}},\mspace{20mu}{v_{{ij}\; 1} = {v_{ij} - v_{ij}^{LB}}},\mspace{20mu}{u_{i\; 1}^{UB} = {u_{i}^{UB} - u_{i}^{LB}}},\mspace{20mu}{v_{{ij}\; 1}^{UB} = {v_{ij}^{UB} - v_{ij}^{LB}}},}} & (67.2) \\ {\mspace{20mu}{{b_{x} = {{\sum\limits_{j = 1}^{N}{R_{ij} \cdot y_{ij}^{LB}}} - {{\overset{\sim}{a}}_{i}{\overset{\sim}{x}}_{i}^{LB}}}},}} & (67.3) \\ {b_{v} = {{R_{ij} \cdot \left( {E_{{ij}\; 1}^{0} + {A_{ij}^{0} \cdot u_{i}^{LB}} + {B_{ij}^{0} \cdot u_{j}^{LB}} + {C_{ij}^{0} \cdot v_{ij}^{LB}} + {D_{ij}^{0} \cdot v_{ji}^{LB}}} \right)} - y_{ij}^{LB}}} & (67.4) \end{matrix}$

According to an embodiment of the invention, the coefficients of the objective function ƒ can be changed: S _(i1) ^(X) =S _(i) ^(X) +V _(i) ^(X) {tilde over (x)} _(i) ^(LB) S _(i1) ^(U) =S _(i) ^(U) +V _(i) ^(U) u _(j) ^(LB) T _(ij1) ^(Y) =T _(ij) ^(Y) +W _(ij) ^(X) {tilde over (x)} _(i) ^(LB) T _(ij1) ^(V) =T _(ij) ^(V) +W _(ij) ^(V) u _(j) ^(LB)  (68) Let v(y _(ij) ,u _(i) ,v _(ij))_(ij) =y _(ij1) −R _(ij)·(A _(ij) ⁰ ·u _(i1) +B _(ij) ⁰ ·u _(j1) +C _(ij) ⁰ ·v _(ij1) +D _(ij) ⁰ ·v _(ji1))−b _(v).  (69) and

$\begin{matrix} {{o\left( {y_{ij},\overset{\sim}{x}} \right)}_{i} = {{{\overset{\sim}{a}}_{i}{\overset{\sim}{x}}_{1}} - {\sum\limits_{j = 1}^{N}{R_{ij} \cdot y_{{ij}\; 1}}} - {b_{x}.}}} & (70) \end{matrix}$

According to an embodiment of the invention, notation can be simplified and readability can be improved by redefining the variable vector x and slack variables s so that

$\begin{matrix} {x = {{\begin{bmatrix} \overset{\sim}{x} \\ u \\ y \\ v \end{bmatrix}\mspace{14mu}{and}\mspace{14mu} s} = \begin{bmatrix} s_{\overset{\sim}{x}} \\ s_{u} \\ s_{y} \\ s_{v} \end{bmatrix}}} & (71) \end{matrix}$

By introducing logarithmic barrier terms, the following Lagrangian function can be derived: L _(μ)(x,s)=ƒ(x)−μ^(T)(ln x+ln s)−η^(T)(x+s−x ^(UB))−λ₁ ^(T) o(y _(ij) ,{tilde over (x)})−λ₂ ^(T) v(y _(ij) ,u _(i) ,v _(ij)).  (72) The Lagrange multiplier vectors η, λ₁, λ₂, are referred to as dual variables.

The KKT optimality conditions of the barrier system are:

$\begin{matrix} {\begin{bmatrix} {\Delta\; L_{X}} \\ {\Delta\; L_{S}} \\ {\Delta\; L_{\eta}} \\ {\Delta\; L_{\lambda_{1}}} \\ {\Delta\; L_{\lambda_{2}}} \end{bmatrix} = {\begin{bmatrix} \begin{matrix} {{\Delta\;{f(x)}} - {\mu/x} - \eta - {{\nabla_{0}\left( {y_{{ij}\; 1},{\overset{\sim}{x}}_{1}} \right)^{T}}\lambda_{1}} -} \\ {v\left( {y_{ij},u_{i},v_{ij}} \right)^{T}\lambda_{2}} \end{matrix} \\ {{{- \mu}/s} - \eta} \\ {- \left( {x + s - x^{UB}} \right)} \\ {- {o\left( {y_{ij},{\overset{\sim}{x}}_{1}} \right)}} \\ {- {v\left( {y_{ij},u_{i},v_{ij}} \right)}} \end{bmatrix} = 0.}} & (73) \end{matrix}$ By defining z_(x)=μ/x and z_(s)=μ/s, EQ. (73) becomes equivalent to:

$\begin{matrix} {\begin{bmatrix} {{\nabla{f(x)}} - z_{x} - \eta - {{\nabla_{0}\left( {y_{{ij}\; 1},{\overset{\sim}{x}}_{1}} \right)^{T}}\lambda_{1}} - {{\nabla{v\left( {y_{ij},u_{i},v_{ij}} \right)}^{T}}\lambda_{2}}} \\ {{- \eta} - z_{s}} \\ {- \left( {x + s - x^{UB}} \right)} \\ {- {o\left( {y_{ij},{\overset{\sim}{x}}_{1}} \right)}} \\ {- {v\left( {y_{ij},u_{i},v_{ij}} \right)}} \\ {{z_{x}^{T}x} - {\mu\; e}} \\ {{z_{s}^{T}s} - {\mu\; e}} \end{bmatrix} = 0} & (74) \end{matrix}$ Newton's method can be used to solve the above non-linear equations, which yields the following linear system of equations, where His a Jacobian matrix:

$\begin{matrix} {{H \cdot \begin{bmatrix} {\Delta\; x} \\ {\Delta\; s} \\ {\Delta\;\eta} \\ {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \\ {\Delta\; z_{x}} \\ {\Delta\; z_{s}} \end{bmatrix}} = {- \begin{bmatrix} \begin{matrix} {{\nabla{f(x)}} - z_{x} - \eta - {{\nabla_{0}\left( {y_{{ij}\; 1},{\overset{\sim}{x}}_{1}} \right)^{T}}\lambda_{1}} -} \\ {{\nabla v}\left( {y_{ij},u_{i},v_{ij}} \right)^{T}\lambda_{2}} \end{matrix} \\ {{- \eta} - z_{s}} \\ {- \left( {x + s - x^{UB}} \right)} \\ {- {o\left( {y_{ij},{\overset{\sim}{x}}_{1}} \right)}} \\ {- {v\left( {y_{ij},u_{i},v_{ij}} \right)}} \\ {{z_{x}^{T}x} - {\mu\; e}} \\ {{z_{s}^{T}s} - {\mu\; e}} \end{bmatrix}}} & (75) \end{matrix}$

Let Z_(x), Z_(s), X and S be the diagonal matrices of z_(x), z_(s), x and s, respectively. Then, the Jacobian matrix H takes the following form:

$\begin{matrix} {\mspace{20mu}{{H = \begin{bmatrix} H_{11} & 0 & {- I} & H_{41}^{T} & H_{51}^{T} & {- I} & 0 \\ 0 & 0 & {- I} & 0 & 0 & 0 & {- I} \\ {- I} & {- I} & 0 & 0 & 0 & 0 & 0 \\ H_{41} & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{51} & 0 & 0 & 0 & 0 & 0 & 0 \\ Z_{x} & 0 & 0 & 0 & 0 & X & 0 \\ 0 & Z_{s} & 0 & 0 & 0 & 0 & S \end{bmatrix}}\mspace{20mu}{where}}} & (76) \\ {\mspace{20mu}{{H_{11} = {{diag}\left( {V^{X},V^{U},W^{Y},W^{V}} \right)}},}} & (77.1) \\ {\mspace{20mu}{{H_{41} = \left\lbrack {{- I},0,R,0} \right\rbrack},}} & (77.2) \\ {\mspace{20mu}{{H_{51} = \left\lbrack {0,H_{51}^{u},{- I},H_{51}^{v}} \right\rbrack},}} & (77.3) \\ {\mspace{20mu}{{R = \begin{pmatrix} \left\lbrack {R_{11}\mspace{14mu}\ldots\mspace{14mu} R_{1N}} \right\rbrack & 0 & \ldots & 0 \\ 0 & \left\lbrack {R_{21}\mspace{14mu}\ldots\mspace{14mu} R_{2N}} \right\rbrack & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & \left\lbrack {R_{N\; 1}\mspace{14mu}\ldots\mspace{14mu} R_{NN}} \right\rbrack \end{pmatrix}},}} & (77.4) \\ {\left( H_{51}^{u} \right)_{ij} = \left\{ \begin{matrix} \begin{pmatrix} 0 & \ldots & {R_{ij}A_{ij}^{0}} & 0 & \ldots & {R_{ij}B_{ij}^{0}} & 0 & \ldots \end{pmatrix} & {{{if}\mspace{14mu} i} < j} \\ \begin{pmatrix} 0 & \ldots & {R_{ij}B_{ij}^{0}} & 0 & \ldots & {R_{ij}A_{ij}^{0}} & 0 & \ldots \end{pmatrix} & {{{if}\mspace{14mu} i} > j} \\ \begin{pmatrix} 0 & \ldots & {{R_{ij}A_{ij}^{0}} + {R_{ij}B_{ij}^{0}}} & 0 & \ldots & 0 & 0 & \ldots \end{pmatrix} & {{{if}\mspace{14mu} i} = j} \end{matrix} \right.} & (77.5) \\ {\left( H_{51}^{v} \right)_{ij} = \left\{ \begin{matrix} \begin{pmatrix} 0 & \ldots & {R_{ij}C_{ij}^{0}} & 0 & \ldots & {R_{ij}D_{ij}^{0}} & 0 & \ldots \end{pmatrix} & {{{if}\mspace{14mu} i} < j} \\ \begin{pmatrix} 0 & \ldots & {R_{ij}D_{ij}^{0}} & 0 & \ldots & {R_{ij}C_{ij}^{0}} & 0 & \ldots \end{pmatrix} & {{{if}\mspace{14mu} i} > j} \\ \begin{pmatrix} 0 & \ldots & {{R_{ij}C_{ij}^{0}} + {R_{ij}D_{ij}^{0}}} & 0 & \ldots & 0 & 0 & \ldots \end{pmatrix} & {{{if}\mspace{14mu} i} = j} \end{matrix} \right.} & \; \end{matrix}$ where (H₅₁ ^(u))_(ij) or (H₅₁ ^(v))_(ij) is the row corresponding to branch (i,j). Note that all four diagonal blocks in the matrix H₁₁ are themselves diagonal by design.

According to an embodiment of the invention, one can define

$\begin{matrix} {{\begin{bmatrix} {{\nabla{f(x)}} - z_{x} - \eta - {{\nabla_{0}\left( {y_{{ij}\; 1},{\overset{\sim}{x}}_{1}} \right)^{T}}\lambda_{1}} - {{\nabla{v\left( {y_{ij},u_{i},v_{ij}} \right)}^{T}}\lambda_{2}}} \\ {{- \eta} - z_{s}} \\ {- \left( {x + s - x^{UB}} \right)} \\ {- {o\left( {y_{ij},{\overset{\sim}{x}}_{1}} \right)}} \\ {- {v\left( {y_{ij},u_{i},v_{ij}} \right)}} \\ {{z_{x}^{T}x} - {\mu\; e}} \\ {{z_{s}^{T}s} - {\mu\; e}} \end{bmatrix} = \begin{bmatrix} r_{x} \\ r_{s} \\ r_{\eta} \\ r_{\lambda_{1}} \\ r_{\lambda_{2}} \\ r_{z_{x}} \\ r_{z_{s}} \end{bmatrix}}\mspace{20mu}{{Then},}} & (78) \\ {\mspace{20mu}{{\begin{bmatrix} H_{11} & 0 & {- I} & H_{41}^{T} & H_{51}^{T} & {- I} & 0 \\ 0 & 0 & {- I} & 0 & 0 & 0 & {- I} \\ {- I} & {- I} & 0 & 0 & 0 & 0 & 0 \\ H_{41} & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{51} & 0 & 0 & 0 & 0 & 0 & 0 \\ Z_{x} & 0 & 0 & 0 & 0 & X & 0 \\ 0 & Z_{s} & 0 & 0 & 0 & 0 & S \end{bmatrix}\begin{bmatrix} {\Delta\; x} \\ {\Delta\; s} \\ {\Delta\;\eta} \\ {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \\ {\Delta\; z_{x}} \\ {\Delta\; z_{s}} \end{bmatrix}} = {- {\begin{bmatrix} r_{x} \\ r_{s} \\ r_{\eta} \\ r_{\lambda_{1}} \\ r_{\lambda_{2}} \\ r_{z_{x}} \\ r_{z_{s}} \end{bmatrix}.}}}} & (79) \end{matrix}$ EQS. (79) are analogous to EQS. (26), above, and like EQS. (26), the dimension of the system of EQS. (79) can be very large. However, a procedure according to an embodiment of the invention that uses a series of back substitutions, described below, similar to that used for solving EQS. (26), can reduce the dimensionality of the system so that the equations require much less computational time to solve.

According to an embodiment of the invention, the above system of Newton equations can be solved by reducing the dimension as follows. From the last two equation there is: Δz _(X) =X ⁻¹(−r _(Z) _(x) −Z _(X) Δx)  (80) Δz _(S) =S ⁻¹(−r _(Z) _(s) −Z _(S) Δs)  (81) EQS. (80) and (81) can be substituted back into the other equations, which can be re-organized into primal and dual parts as follows:

$\begin{matrix} {{\begin{bmatrix} H_{41} & 0 & | & 0 & 0 & 0 \\ H_{51} & 0 & | & 0 & 0 & 0 \\ {- I} & {- I} & | & 0 & 0 & 0 \\  - & - & \; & - & - & - \\ {H_{11} + {X^{- 1}Z_{x}}} & 0 & | & {- I} & H_{41}^{T} & H_{51}^{T} \\ 0 & {S^{- 1}Z_{s}} & | & {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; x} \\ {\Delta\; s} \\  - \\ {\Delta\;\eta} \\ {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \end{bmatrix}} = \begin{bmatrix} {- r_{\lambda_{1}}} \\ {- r_{\lambda_{2}}} \\ {- r_{\eta}} \\  - \\ {{- r_{x}} - {X^{- 1}r_{z_{x}}}} \\ {{- r_{s}} - {S^{- 1}r_{z_{s}}}} \end{bmatrix}} & (82) \end{matrix}$

According to an embodiment of the invention, define:

$\begin{matrix} {{{\overset{\sim}{H}}_{11} = {H_{11} + {X^{- 1}Z_{x}}}}{{\overset{\sim}{H}}_{22} = {S^{- 1}Z_{s}}}} & (83.1) \\ {{{\overset{\sim}{r}}_{x} = {r_{x} + {X^{- 1}r_{z_{x}}}}}{{\overset{\sim}{r}}_{s} = {r_{s} + {S^{- 1}r_{Z_{s}}}}}} & (83.2) \\ {{\Delta\; x} = \begin{bmatrix} \left\{ {\Delta\;\overset{\sim}{x}} \right\} \\ \left\{ {\Delta\; y} \right\} \\ \left\{ {\Delta\; u} \right\} \\ \left\{ {\Delta\; v} \right\} \end{bmatrix}} & (83.3) \\ {{{\Delta\;\lambda_{1}} = \left\{ {\Delta\;\lambda_{1}^{\overset{\sim}{x}}} \right\}}{{\Delta\;\lambda_{2}} = \left\{ {\Delta\;\lambda_{2}^{y}} \right\}}} & (83.4) \end{matrix}$ Therefore, according to an embodiment of the invention, one has:

$\begin{matrix} {{{\begin{bmatrix} H_{41} & 0 \\ H_{51} & 0 \\ {- I} & {- I} \end{bmatrix}\begin{bmatrix} {\Delta\; x} \\ {\Delta\; s} \end{bmatrix}} = \begin{bmatrix} {- r_{\lambda_{1}}} \\ {- r_{\lambda_{2}}} \\ {- r_{\eta}} \end{bmatrix}},{and}} & (84) \\ {{{\begin{bmatrix} {\overset{\sim}{H}}_{11} & 0 \\ 0 & {\overset{\sim}{H}}_{22} \end{bmatrix}\begin{bmatrix} {\Delta\; x} \\ {\Delta\; s} \end{bmatrix}} + {\begin{bmatrix} {- I} & H_{41}^{T} & H_{51}^{T} \\ {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\;\eta} \\ {\Delta\;\lambda_{1}} \\ {\Delta\;\lambda_{2}} \end{bmatrix}}} = {\begin{bmatrix} {- {\overset{\sim}{r}}_{x}} \\ {- {\overset{\sim}{r}}_{s}} \end{bmatrix}.}} & (85) \end{matrix}$ Back Substitution Solution Series for Reduced Newton System

More specifically, similar to the solution of EQS. (30) and (31), above, EQS. (84) and (85) can be expanded using the definitions of EQS. (83.1)-(83.4) and the model definitions of EQS. (67.1)-(67.4), (68)-(71), and (77.1)-(77.5). Thus, EQ. (84) has the following expanded form for each bus i and branch (i, j):

$\begin{matrix} {\mspace{20mu}{{{{{\overset{\sim}{a}}_{i}\Delta\;{\overset{\sim}{x}}_{1}} - {\sum\limits_{j}{R_{ij}\Delta\; y_{{ij}\; 1}}}} = \left( r_{\lambda_{1}} \right)_{i}}{{{\Delta\; y_{{ij}\; 1}} - {R_{ij}\left( {{A_{ij}^{0}\Delta\; u_{i\; 1}} + {B_{ij}^{0}\Delta\; u_{j\; 1}} + {C_{ij}^{0}\Delta\; v_{{ij}\; 1}} + {D_{ij}^{0}\Delta\; v_{{ji}\; 1}}} \right)}} = \left( r_{\lambda_{2\;}} \right)_{ij}}\mspace{20mu}{{{\Delta\; x} + {\Delta\; s}} = r_{\eta}}}} & (86) \end{matrix}$ Similarly, EQ. (85) has the following expanded form:

$\begin{matrix} {\mspace{20mu}{{{{\Delta\;\eta} - {\Delta\;\lambda_{1}}} = {- {\overset{\sim}{r}}_{\beta\; 1}}}\mspace{20mu}{{{{\overset{\sim}{V}}^{X}\Delta\;{\overset{\sim}{x}}_{1}} - {\Delta\;\eta} - {{\overset{\sim}{A}}^{T}\Delta\;\lambda_{1}}} = {- {\overset{\sim}{r}}_{{\overset{\sim}{x}}_{1}}}}}} & (87.1) \\ {\mspace{20mu}{{{{{\overset{\sim}{W}}^{Y}\Delta\; y_{1}} - {\Delta\;\eta} + {R\;\Delta\;\lambda_{1}} - {\Delta\;\lambda_{2}}} = {- {\overset{\sim}{r}}_{y_{1}}}}\mspace{20mu}{{{{\overset{\sim}{H}}_{22}\Delta\; s} - {\Delta\;\eta}} = {- {\overset{\sim}{r}}_{s}}}\mspace{20mu}{and}}} & (87.2) \\ {{{{{\overset{\sim}{V}}_{i}^{U}\Delta\; u_{i\; 1}} - {\Delta\;\eta_{i}} + {\sum\limits_{j}{R_{ij}A_{ij}^{T}\Delta\;\lambda_{2}^{y_{{ij}\; 1}}}} + {\sum\limits_{j}{R_{ji}B_{ji}^{T}\Delta\;\lambda_{2}^{y_{{ji}\; 1}}}}} = {- \left( {\overset{\sim}{r}}_{u_{1}} \right)_{i}}}\mspace{20mu}{{{{\overset{\sim}{W}}_{ij}^{V}\Delta\; v_{{ij}\; 1}} - {\Delta\;\eta_{ij}}\; + {R_{ij}C_{ij}^{T}\Delta\;\lambda_{2}^{y_{{ij}\; 1}}} + {R_{ji}D_{ji}^{T}\Delta\;\lambda_{2}^{v_{{ji}\; 1}}}} = {- \left( {\overset{\sim}{r}}_{v_{1}} \right)_{ij}}}} & (87.3) \end{matrix}$ Suppose there is a candidate Newton step solution for the control variables ΔU=[{Δu_(i1)}, {Δv_(ij1)}]^(T). Then, one can solve for all the other variables directly through EQS. (86) and (87.1)-(87.2) forming a back substitution series of direct linear elimination: {Δu ₁ },{Δv ₁ }→{Δy ₁}→{Δ{tilde over (x)}₁ }→Δs→Δη→Δλ ₁→Δλ₂  (88) Then, this back substitution series rule of EQ. (88) can be used to define a linear operator L:ΔU→{Δy₁, Δ{tilde over (x)}₁, ΔS, Δη, Δλ₁, Δλ₂} defined over the span of ΔU=[{Δu_(i1)}, {Δv_(ij1)}]^(T). This operator can be substituted back into the linear system of EQS. (87.3) which then becomes a system of joint linear equations for the ΔU only. This system is typically full-rank, well-defined and its size is a few orders of magnitude smaller than the reduced Newton system of EQS. (85). In the simplest case of {Δv_(ij1)}=∅, this system is N by N and relatively sparse (density factor ˜0.1%) since the non-zero entries for any row i are given by the set of “double-neighbor” nodes S₂(i)={j: ∃k s.t. R_(jk)=1, R_(k1)=1}. This system is best solved via indirect or iterative methods, such as the Preconditioned Conjugate Method, since computing the full matrix form of the linear operator L is too costly for massive-scale problems. System Implementation

It is to be understood that embodiments of the present invention can be implemented in various forms of hardware, software, firmware, special purpose processes, or a combination thereof. In one embodiment, the present invention can be implemented in software as an application program tangible embodied on a computer readable program storage device. The application program can be uploaded to, and executed by, a machine comprising any suitable architecture.

FIG. 4 is a block diagram of an exemplary computer system for implementing an interior point method solution for an optimal power flow model of a smart power grid according to an embodiment of the invention. Referring now to FIG. 4, a computer system 41 for implementing the present invention can comprise, inter alia, a central processing unit (CPU) 42, a memory 43 and an input/output (I/O) interface 44. The computer system 41 is generally coupled through the I/O interface 44 to a display 45 and various input devices 46 such as a mouse and a keyboard. The support circuits can include circuits such as cache, power supplies, clock circuits, and a communication bus. The memory 43 can include random access memory (RAM), read only memory (ROM), disk drive, tape drive, etc., or a combinations thereof. The present invention can be implemented as a routine 47 that is stored in memory 43 and executed by the CPU 42 to process the signal from the signal source 48. As such, the computer system 41 is a general purpose computer system that becomes a specific purpose computer system when executing the routine 47 of the present invention.

The computer system 41 also includes an operating system and micro instruction code. The various processes and functions described herein can either be part of the micro instruction code or part of the application program (or combination thereof) which is executed via the operating system. In addition, various other peripheral devices can be connected to the computer platform such as an additional data storage device and a printing device.

It is to be further understood that, because some of the constituent system components and method steps depicted in the accompanying figures can be implemented in software, the actual connections between the systems components (or the process steps) may differ depending upon the manner in which the present invention is programmed. Given the teachings of the present disclosure provided herein, one of ordinary skill in the related art will be able to contemplate these and similar implementations or configurations of the present invention.

While embodiments of the present invention has been described in detail with reference to exemplary embodiments, those skilled in the art will appreciate that various modifications and substitutions can be made thereto without departing from the spirit and scope of the disclosure as set forth in the appended claims. 

What is claimed is:
 1. A method for approximating an optimal power flow of a smart electric power grid, comprising the steps of: providing a cost function that models a smart electric power grid having N buses iε[1,N] connected by branches i,jε[1,N]×[1, N] with connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\mspace{14mu}{{branch}\mspace{14mu} i}},j,} \\ {0,} & {{otherwise},} \end{matrix} \right.$ in terms of bus control vector u_(i) that include voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(1•), bus dependent variables vector x_(i) that include net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, wherein superscripts P and Q respectively refer to the net active power injection and the reactive power injection, branch control vector u_(ij), that includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij), that includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, constraining the bus and branch control vectors, and the bus and branch dependent variables vectors based on physical and operational requirements and actual measurements; deriving a set of linear equations that minimize said cost function subject to the constraints, from an expression of an extremum of said cost function with respect to all arguments; reducing a dimension of the linear equations by solving for a subset of the linear equations; re-organizing the reduced dimension linear equations into primal and dual parts; and decomposing the re-organized reduced dimensional linear equations into two systems of block matrix equations, wherein a solution of said two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through said smart electric power grid, wherein the steps of reducing a dimension of the linear equations, re-organizing the reduced dimension linear equations, and decomposing the re-organized reduced dimensional linear equations into two systems of block matrix equations are performed by a computer processor.
 2. The method of claim 1, wherein said cost function is ${{f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\;\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}u_{i}^{T}V_{i}^{U}u_{i}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}u_{i}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}{R_{ij} \cdot \;\left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}u_{ij}^{T}W_{ij}^{U}u_{ij}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}u_{ij}}} \right\rbrack}}}},$ wherein bus control vector u_(i) includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(1•), bus dependent variables vector x_(i) includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active X_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, and connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ wherein ${x_{i} = {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{ij}}}},$ wherein V^(X) is a Hessian of the objective function with respect to x_(i), V^(U) is a Hessian of the objective function with respect to u_(i), S^(X) is a gradient of the objective function with respect to x_(i), S^(U) is a gradient of the objective function with respect to u_(i),W^(X) is a Hessian of the objective function with respect to x_(ij), W^(U) is a Hessian of the objective function with respect to u_(ij), T^(X) is a gradient of the objective function with respect to x_(ij), and T^(U) is a gradient of the objective function with respect to u_(ij), subject to the constraints x _(i) ^(LB) ≦x _(i) ≦x _(i) ^(UB) , x _(ij) ^(LB) ≦x _(ij) ≦x _(ij) ^(UB), u _(i) ^(LB) ≦u _(i) ≦u _(i) ^(UB) , u _(ij) ^(LB) ≦u _(ij) ≦u _(ij) ^(UB), wherein superscripts LB and UB respectively refer to a lower bound and an upper bound, and satisfying a voltage law x_(ij)=R_(ij)·(A_(ij) ⁰u_(i)+B_(ij) ⁰u_(j)+C_(ij) ⁰u_(ij)+D_(ij) ⁰u_(ji)+E_(ij) ⁰), wherein the matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰,u_(j) ⁰,u_(ij) ⁰).
 3. The method of claim 2, wherein said linear equations that minimize said cost function are: ${{\begin{bmatrix} H_{11} & 0 & {- I} & H_{41}^{T} & H_{51}^{T} & {- I} & 0 \\ 0 & 0 & {- I} & 0 & 0 & 0 & {- I} \\ {- I} & {- I} & 0 & 0 & 0 & 0 & 0 \\ H_{41} & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{51} & 0 & 0 & 0 & 0 & 0 & 0 \\ {{diag}\left( Z_{x} \right)} & 0 & 0 & 0 & 0 & {{diag}(X)} & 0 \\ 0 & {{diag}\left( Z_{s} \right)} & 0 & 0 & 0 & 0 & {{diag}(S)} \end{bmatrix} \cdot \begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; Z_{x}} \\ {\Delta\; Z_{s}} \end{bmatrix}} = {- \begin{bmatrix} {{\nabla\;{f(X)}} - Z_{x} - Y - {y_{1}^{T}{\nabla\;{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla\;{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- \;{v(X)}} \\ {{Z_{x} \cdot X} - {\mu\; e}} \\ {{Z_{s} \cdot S} - {\mu\; e}} \end{bmatrix}}},$ wherein [ΔX, ΔS, ΔY, Δy₁, Δy₂, ΔZ_(x), ΔZ_(s)] are unknowns, I represents an identity matrix, $\mspace{79mu}{{H_{11} = {{diag}\left( {V^{x},W^{x},V^{u},W^{U}} \right)}},\mspace{79mu}{H_{41} = \left\lbrack {{- I},R,0,0} \right\rbrack},{R = \begin{pmatrix} \left\lbrack {R_{11}I\mspace{14mu}\ldots\mspace{14mu} R_{1N\; 1}I} \right\rbrack & 0 & \ldots & 0 \\ 0 & \left\lbrack {R_{21}I\mspace{14mu}\ldots\mspace{14mu} R_{2N\; 2}I} \right\rbrack & \ldots & 0 \\ 0 & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & \left\lbrack {R_{N_{1}1}I\mspace{14mu}\ldots\mspace{14mu} R_{N_{1}N_{1}}I} \right\rbrack \end{pmatrix}},}$ a row of branch (i, j) of H₅₁ is ${\left( H_{51} \right)_{ij} = \left\lbrack {0,{- I},{\ldots\mspace{14mu} R_{ij}A_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}B_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}C_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}D_{ij}^{0}},{\ldots\mspace{14mu} 0}} \right\rbrack},{X = \left( {\left\{ x_{il} \right\},\left\{ x_{ijl} \right\},\left\{ u_{il} \right\},\left\{ u_{ijl} \right\}} \right)^{T}},{{{slack}\mspace{14mu}{variables}\mspace{14mu} S} = \left( {\left\{ s_{xi} \right\},\left\{ s_{xij} \right\},\left\{ s_{ui} \right\},\left\{ s_{uij} \right\}} \right)^{T}},{{o\left( {x_{ij},x_{i}} \right)} = {x_{i\; 1} - {\sum\limits_{j = 1}^{N}\;{R_{ij} \cdot x_{{ij}\; 1}}} - b_{x}}},{{\upsilon(X)} = {x_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot u_{{ij}\; 1}} + {D_{ij}^{0} \cdot u_{{ji}\; 1}}} \right)} - b_{\upsilon}}},$ Y, y₁, and y₂ are dual variables, and Z_(x)=μ/X, Z_(s)=μ/S wherein μ is a barrier parameter.
 4. The method of claim 3, wherein reducing a dimension of the linear equations further comprises solving for ΔZ _(X)=diag(X)⁻¹(−r _(Z) _(X) −diag(Z _(X))ΔX), ΔZ _(S)=diag(S)⁻¹(−r _(Z) _(S) −diag(Z _(S))ΔS), wherein r_(Z) _(X) =Z_(X)·X−μe, r_(Z) _(S) =Z_(S)·S−μe.
 5. The method of claim 4, wherein the reduced dimension linear equations re-organized into primal and dual parts are ${\begin{bmatrix} H_{41} & 0 & ❘ & 0 & 0 & 0 \\ H_{51} & 0 & ❘ & 0 & 0 & 0 \\ {- I} & {- I} & ❘ & 0 & 0 & 0 \\  - & - & ❘ & - & - & - \\ {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{x} \right)}}} & 0 & ❘ & {- I} & H_{41}^{T} & H_{51}^{T} \\ 0 & {{{diag}(S)}^{- 1}{{diag}\left( Z_{s} \right)}} & ❘ & {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\  - \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix}} = {\quad{\begin{bmatrix} {- r_{y\; 1}} \\ {- r_{y\; 2}} \\ {- r_{Y}} \\  - \\ {{- r_{X}} - {{{diag}(X)}^{- 1}{rZ}_{x}}} \\ {{- r_{S}} - {{{diag}(S)}^{- 1}{rZ}_{s}}} \end{bmatrix},\mspace{79mu}{{{wherein}\mspace{79mu}\begin{bmatrix} r_{X} \\ r_{S} \\ r_{Y} \\ r_{y_{1}} \\ r_{y_{2}} \end{bmatrix}} = {\begin{bmatrix} {{\nabla\;{f(X)}} - Z_{x} - Y - {y_{1}^{T}{\nabla\;{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- \;{v(X)}} \end{bmatrix}.}}}}$
 6. The method of claim 5, wherein the two systems of block matrix equations are ${{\begin{bmatrix} H_{41} & 0 \\ H_{51} & 0 \\ {- I} & {- I} \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} = \begin{bmatrix} {- r_{y\; 1}} \\ {- r_{y\; 2}} \\ {- r_{Y}} \end{bmatrix}},{{{{{and}\begin{bmatrix} {\overset{\sim}{H}}_{11} & 0 \\ 0 & {\overset{\sim}{H}}_{22} \end{bmatrix}}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} + {\begin{bmatrix} {- I} & H_{41}^{T} & H_{51}^{T} \\ {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix}}} = \begin{bmatrix} {- {\overset{\sim}{r}}_{X}} \\ {- {\overset{\sim}{r}}_{S}} \end{bmatrix}},{wherein}$ ${{\overset{\sim}{H}}_{11} = {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{X} \right)}}}},{{\overset{\sim}{H}}_{22} = {{{diag}(S)}^{- 1}{{diag}\left( Z_{S} \right)}}},{{\overset{\sim}{r}}_{X} = {r_{X} + {{{diag}(X)}^{- 1}r_{Z_{X}}}}},{{\overset{\sim}{r}}_{S} = {r_{S} + {{{diag}(S)}^{- 1}r_{Z_{S}}}}},{{\Delta\; X} = \begin{bmatrix} \left\{ {\Delta\; x_{i}} \right\} \\ \left\{ {\Delta\; x_{ij}} \right\} \\ \left\{ {\Delta\; u_{i}} \right\} \\ \left\{ {\Delta\; u_{ij}} \right\} \end{bmatrix}},{{\Delta\; y_{1}} = \left\{ {\Delta\; y_{1}^{x_{i}}} \right\}},{{\Delta\; y_{2}} = {\left\{ {\Delta\; y_{2}^{x_{ij}}} \right\}.}}$
 7. The method of claim 6, wherein the first block matrix equation is expanded for each bus i and branch(i, j) as follows: ${{\Delta\; x_{i}} - {\sum\limits_{j}^{\;}\;{R_{ij}\Delta\; x_{ij}}}} = r_{y_{1}}^{x_{i}}$ Δ x_(ij) − R_(ij)(A_(ij)Δ u_(i) + B_(ij)Δ u_(j) + C_(ij)Δ u_(ij) + D_(ij)Δ u_(ji)) = r_(y₂)^(x_(ij)) Δ X + Δ S = r_(Y), and the second block matrix equation is expanded for each bus i and branch(i, j) as follows: ${{{\overset{\sim}{V}}_{i}^{X}\Delta\; x_{i}} - {\Delta\; Y^{x_{i}}} - {\Delta\; y_{1}^{x_{i}}}} = {- {\overset{\sim}{r}}_{x_{i}}}$ ${{{\overset{\sim}{W}}_{ij}^{X}\Delta\; x_{ij}} - {\Delta\; Y^{x_{ij}}} + {R_{ij}\Delta\; y_{1}^{x_{i}}} - {\Delta\; y_{2}^{x_{ij}}}} = {- {\overset{\sim}{r}}_{x_{ij}}}$ ${{{\overset{\sim}{V}}_{i}^{U}\Delta\; u_{i}} - {\Delta\; Y^{u_{i}}} + {\sum\limits_{j}^{\;}\;{R_{ij}A_{ij}^{T}\Delta\; y_{2}^{x_{ij}}}} + {\sum\limits_{j}^{\;}\;{R_{ji}B_{ji}^{T}\Delta\; y_{2}^{x_{ji}}}}} = {- {\overset{\sim}{r}}_{u_{i}}}$ ${{{\overset{\sim}{W}}_{ij}^{U}\Delta\; u_{ij}} - {\Delta\; Y^{u_{ij}}} + {R_{ij}C_{ij}^{T}\Delta\; y_{2}^{x_{ij}}} + {R_{ji}D_{ji}^{T}\Delta\; y_{2}^{x_{ji}}}} = {- {\overset{\sim}{r}}_{u_{ij}}}$ ${{{{\overset{\sim}{H}}_{22}\Delta\; S} - {\Delta\; Y}} = {- {\overset{\sim}{r}}_{S}}},$ and further comprising solving for the other variables in the block matrix equations in the following order: {Δu _(i) },{Δu _(ij) }→{Δx _(ij) }→{Δx _(i) }→ΔS→ΔY→Δy ₁ →Δy ₂.
 8. The method of claim 6, further comprising: initializing unknowns [ΔX^(k), ΔS^(k), ΔY^(k), Δy₁ ^(k), Δy₂ ^(k)] wherein k is an iteration counter; updating the barrier parameter as μ^(k+1)←ε·μ^(k), εε(0.1, . . . , 0.9); evaluating (ΔX^(k), ΔS^(k), ΔY^(k), Δy₁ ^(k), Δy₂ ^(k)) using the two block matrix equations; calculating (ΔZ_(X) ^(k),ΔZ_(S) ^(k)); determining a step length Δβ^(k) from a line search; and updating X ^(k+1) ←X ^(k)+Δβ^(k) ΔX ^(k), S ^(k+1) ←S ^(k)+Δβ^(k) ΔS ^(k), Y ^(k+1) ←Y ^(k)+Δβ^(k) ΔY ^(k), y ₁ ^(k+1) ←y ₁ ^(k)+Δβ^(k) Δy ₁ ^(k), y ₂ ^(k+1) ←y ₂ ^(k)+Δβ^(k) Δy ₂ ^(k), Z _(X) ^(k+1) ←Z _(X) ^(k)+Δβ^(k) ΔZ _(X) ^(k), Z _(S) ^(k+1) ←Z _(S) ^(k)+Δβ^(k) ΔZ _(S) ^(k), wherein the steps of updating the barrier parameter, evaluating (ΔX^(k), ΔS^(k), ΔY^(k), Δy₁ ^(k), Δy₂ ^(k)), calculating (ΔZ_(X) ^(k),ΔZ_(S) ^(k)), determining a step length Δβ^(k), and updating X^(k+1), S^(k+1), Y^(k+1), y₁ ^(k+1), y₂ ^(k+1), Z_(X) ^(k+1) and Z_(S) ^(k+1) are repeated until either a dual gap or a maximum of (r_(X), r_(S)) is less than a pre-specified minimum.
 9. The method of claim 1, wherein said cost function is, ${f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\;\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}{\sum\limits_{Ki}{u_{i}^{kT}V_{i}^{U}u_{i}^{k}}}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}{\sum\limits_{Ki}u_{i}^{k}}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}{R_{ij} \cdot \left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}{\sum\limits_{Kij}{u_{ij}^{kT}W_{ij}^{U}u_{ij}^{k}}}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}{\sum\limits_{Kik}u_{ij}^{k}}}} \right\rbrack}}}$ wherein bus control vector u_(i) includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, and connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ wherein ${x_{i} = {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{ij}}}},$ wherein V^(X) is a Hessian of the objective function with respect to x_(i), V^(U) is a Hessian of the objective function with respect to u_(i), S^(X) is a gradient of the objective function with respect to x_(i), S^(U) is a gradient of the objective function with respect to u_(i), W^(X) is a Hessian of the objective function with respect to x_(ij), W^(U) is a Hessian of the objective function with respect to u_(ij), T^(X) is a gradient of the objective function with respect to x_(ij), and T^(U) is a gradient of the objective function with respect to u_(ij), subject to constraints ${x_{i}^{LB} \leq x_{i} \leq x_{i}^{UB}},{u_{i}^{LB} \leq {\sum\limits_{Ki}u_{i}^{k}} \leq u_{i}^{UB}},{x_{ij}^{LB} \leq x_{ij} \leq x_{ij}^{UB}},{u_{ij}^{LB} \leq {\sum\limits_{Kij}u_{ij}^{k}} \leq u_{ij}^{UB}},$ wherein a control grid is discretized by binary variables zε{0,1} as z _(i) ^(k) u _(ik) ^(LB) ≦u _(i) ^(k) ≦z _(i) ^(k) u _(ik) ^(UB), Z _(ij) ^(k) u _(ijk) ^(LB) ≦u _(ij) ^(k) ≦z _(ij) ^(k) u _(ijk) ^(UB), wherein the binary variables are subject to constraints $z_{i}^{LB} \leq {\sum\limits_{Ki}z_{i}^{k}} \leq z_{i}^{UB}$ $z_{i}^{LB} \leq {\sum\limits_{Ki}z_{i}^{k}} \leq z_{i}^{UB}$ and satisfying a voltage law x_(ij)=R_(ij)·(Σ_(Ki)A_(ij) ^(k)u_(i) ^(k)+Σ_(Kj)B_(ij) ^(k)u_(j) ^(k)+Σ_(Kij)C_(ij) ^(k)u_(ij) ^(k)+Σ_(Kji)D_(ij) ^(k)u_(ji) ^(k)), where the matrices A_(ij) ⁰,B_(ij) ⁰,C_(ij) ⁰,D_(ij) ⁰,E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰,u_(j) ⁰,u_(ij) ⁰).
 10. The method of claim 1, wherein said cost function is ƒ({tilde over (x)},u,y,v)=½[{tilde over (x)}^(T)V^(X){tilde over (x)}+u^(T)V^(U)u+y^(T)W^(Y)y+v^(T)W^(V)v]+S^(X){tilde over (x)}+S^(U)u+T^(Y)y+T^(V)v, wherein bus control vector u includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector v includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector y includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, wherein bus dependent vector is defined as ${x = {\left. {\overset{\sim}{A}\;\overset{\sim}{x}}\Leftrightarrow x_{i} \right. = {{\sum\limits_{k = 1}^{G_{i}}x_{i,k}^{g}} - {\sum\limits_{k = 1}^{L_{i}}x_{i,k}^{l}} + {\sum\limits_{k = 1}^{{SH}_{i}}x_{i,k}^{sh}}}}},$ wherein where G_(i), L_(i) and SH_(i) are the number of generators, loads and shunts placed on bus i and x_(i,k) ^(g), x_(i,k) ^(i) and x_(i,k) ^(sh) are the k^(th) active and reactive generation, load and shunt variables at bus i, respectively, and Ãε{−1,0,1}^(N×2(G+L+SH)), wherein V^(X) is a Hessian of the objective function with respect to x, V^(U) is a Hessian of the objective function with respect to u, S^(X) is a gradient of the objective function with respect to x, S^(U) is a gradient of the objective function with respect to u, W^(Y) is a Hessian of the objective function with respect to y, W^(V) is a Hessian of the objective function with respect to v, T^(Y) is a gradient of the objective function with respect to y, and T^(V) is a gradient of the objective function with respect to v, said cost function having connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ wherein ${a_{i}\overset{\sim}{x}} = {\sum\limits_{j = 1}^{N}{R_{ij}y_{ij}}}$ wherein a_(i) is i^(th) row of matrix Ã, subject to the constraints {tilde over (x)} ^(LB) ≦{tilde over (x)}≦{tilde over (x)} ^(UB) ,y ^(LB) ≦y≦y ^(UB), u ^(LB) ≦u≦u ^(UB) ,v ^(LB) ≦v≦v ^(UB), and satisfying a voltage law y_(ij)=R_(ij)·(A_(ij) ⁰·u_(i)+β_(ij) ⁰·u_(j)+C_(ij) ⁰·v_(ij)+D_(ij) ⁰·v_(ji)+E_(ij) ⁰), wherein the matrices A_(ij) ⁰,B_(ij) ⁰,C_(ij) ⁰,D_(ij) ⁰,E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰,u_(j) ⁰,u_(ij) ⁰).
 11. A method for approximating an optimal power flow of a smart electric power grid, comprising the steps of: providing a cost function that models a smart electric power grid having N buses iε[1,N] connected by branches i,jε[1,N]×[1,N] with connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\mspace{14mu}{{branch}\mspace{14mu} i}},j,} \\ {0,} & {{otherwise},} \end{matrix} \right.$ in terms of bus control vector u_(i) that include voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(1•), bus dependent variables vector x_(i) that include net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, wherein superscripts P and Q respectively refer to the net active power injection and the reactive power injection, branch control vector u_(ij), that includes transformer tap selections t_(ij) and phase shifts φ_(i), and dependent branch variable vector x_(ij), that includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, constraining the bus and branch control vectors, and the bus and branch dependent variables vectors based on physical and operational requirements and actual measurements; deriving a set of non-linear vector equations $\overset{\sim}{L} = {\begin{bmatrix} {{\nabla{f(X)}} - {\mu/X} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{{- \mu}/S} - Y} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i\;}} \right)}} \\ {- {v(X)}} \end{bmatrix} = 0}$ that minimize said cost function subject to the constraints, from an expression of an extremum of said cost function with respect to all arguments, wherein superscript T refers to a vector or matrix transpose, wherein X=({x_(i)}, {x_(ij)}, {u_(i)}, {u_(ij)})^(T), bus control vector u_(i) includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, S=({s_(xi)}, {s_(xij)}, {s_(ui)}, {s_(uij)})^(T), where s_(xi), s_(ui), s_(xij), s_(uij) are slack variables, ${{o\left( {x_{ij},x_{i\;}} \right)} = {x_{i\; 1} - {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{{ij}\; 1}}} - b_{x}}},{{\upsilon(X)} = {x_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot u_{{ij}\; 1}} + {D_{ij}^{0} \cdot u_{{ji}\; 1}}} \right)} - b_{\upsilon}}}$ b_(x) and b_(v) are predetermined constants, $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ is a connection matrix wherein ${x_{i} = {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{ij}}}},{x_{ij} = {R_{ij} \cdot \left( {{A_{ij}^{0}u_{i}} + {B_{ij}^{0}u_{j}} + {C_{ij}^{0}u_{ij}} + {D_{ij}^{0}u_{ji}} + E_{ij}^{0}} \right)}}$ is a voltage law wherein the matrices A_(ij) ⁰, B_(ij) ⁰, C_(ij) ⁰, D_(ij) ⁰, E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰,u_(j) ⁰,u_(ij) ⁰) Y, y₁, and y₂ are dual variables, and μ is a barrier parameter; re-organizing the vector {tilde over (L)} as ${\overset{\sim}{L} = \begin{bmatrix} \left\{ L^{i} \right\} \\ \left\{ L^{ij} \right\} \end{bmatrix}},$ wherein each L^(i) represent the constraints of bus i, and each L^(ij) represents the constraints of branch (i, j); applying Newton's method to every bus L^(i) and branch L^(ij), to obtain ${{{H^{i} \cdot \begin{bmatrix} {\Delta\; X^{i}} \\ {\Delta\; S^{i}} \\ {\Delta\; Y^{i}} \\ {\Delta\; y_{1}^{i}} \\ {\Delta\; y_{2}^{i}} \end{bmatrix}} = {- L^{i}}};{{H^{ij} \cdot \begin{bmatrix} {\Delta\; X^{ij}} \\ {\Delta\; S^{ij}} \\ {\Delta\; Y^{ij}} \\ {\Delta\; y_{1}^{ij}} \\ {\Delta\; y_{2}^{ij}} \end{bmatrix}} = {- L^{ij}}}},$ where H^(i) and H^(ij) are Jacobian matrices, ΔX^(i)=({Δx_(i)}, {Δu_(i)})^(T), ΔX^(ij)=({Δx_(ij)}, {Δu_(ij)}), ΔS^(i)=({Δs_(xi)}, {Δs_(ui)})^(T), ΔS^(ij)=({Δs_(xij)}, {Δs_(uij)}), and [ΔX, ΔS, ΔY, Δy₁, Δy₂] are unknowns; and using cutting planes to solve these block-wise Newton's equations, wherein a solution of said two systems of block matrix equations yields conditions for a lowest cost per kilowatt hour delivered through said smart electric power grid, wherein the steps of re-organizing the vector {tilde over (L)}, applying Newton's method, and using cutting planes to solve these block-wise Newton's equations, are performed by a computer processor.
 12. A non-transitory program storage device readable by a computer, tangibly embodying a program of instructions executed by the computer to perform the method steps for approximating an optimal power flow of a smart electric power grid, the method comprising the steps of: providing a cost function that models a smart electric power grid having N buses iε[1, N] connected by branches i,jε[1,N]×[1,N] with connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\mspace{14mu}{{branch}\mspace{14mu} i}},j,} \\ {0,} & {{otherwise},} \end{matrix} \right.$ in terms of bus control vector u_(i) that include voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(1•), bus dependent variables vector x_(i) that include net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, wherein superscripts P and Q respectively refer to the net active power injection and the reactive power injection, branch control vector u_(ij), that includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij), that includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, constraining the bus and branch control vectors, and the bus and branch dependent variables vectors based on physical and operational requirements and actual measurements; deriving a set of linear equations that minimize said cost function subject to the constraints, from an expression of an extremum of said cost function with respect to all arguments; reducing a dimension of the linear equations by solving for a subset of the linear equations; re-organizing the reduced dimension linear equations into primal and dual parts; and decomposing the re-organized reduced dimensional linear equations into two systems of block matrix equations, wherein a solution of said two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through said smart electric power grid.
 13. The computer readable program storage device of claim 12, wherein said cost function is ${{f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}u_{i}^{T}V_{i}^{U}u_{i}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}u_{i}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}{R_{ij} \cdot \left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}u_{ij}^{T}W_{ij}^{U}u_{ij}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}u_{ij}}} \right\rbrack}}}},$ wherein bus control vector u_(i) includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, and connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ wherein ${x_{i} = {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{ij}}}},$ wherein V^(X) is a Hessian of the objective function with respect to x_(i), V^(U) is a Hessian of the objective function with respect to u_(i), S^(X) is a gradient of the objective function with respect to x_(i), S^(U) is a gradient of the objective function with respect to u_(i), W^(X) is a Hessian of the objective function with respect to x_(ij), W^(U) is a Hessian of the objective function with respect to u_(ij), T^(X) is a gradient of the objective function with respect to x_(ij), and T^(U) is a gradient of the objective function with respect to u_(ij), subject to the constraints x _(i) ^(LB) ≦x _(i) ≦x _(i) ^(UB) , x _(ij) ^(LB) ≦x _(ij) ≦x _(ij) ^(UB), u _(i) ^(LB) ≦u _(i) ≦u _(i) ^(UB) , u _(ij) ^(LB) ≦u _(ij) ≦u _(ij) ^(UB), wherein superscripts LB and UB respectively refer to a lower bound and an upper bound, and satisfying a voltage law x_(ij)=R_(ij)·(A_(ij) ⁰+B_(ij) ⁰u_(j)+C_(ij) ⁰u_(ij)+D_(ij) ⁰u_(ji)+E_(ij) ⁰), wherein the matrices A_(ij) ⁰,B_(ij) ⁰,C_(ij) ⁰,D_(ij) ⁰,E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰,u_(j) ⁰,u_(ij) ⁰).
 14. The computer readable program storage device of claim 13, wherein said linear equations that minimize said cost function are: ${{\begin{bmatrix} H_{11} & 0 & {- I} & H_{41}^{T} & H_{51}^{T} & {- I} & 0 \\ 0 & 0 & {- I} & 0 & 0 & 0 & {- I} \\ {- I} & {- I} & 0 & 0 & 0 & 0 & 0 \\ H_{41} & 0 & 0 & 0 & 0 & 0 & 0 \\ H_{51} & 0 & 0 & 0 & 0 & 0 & 0 \\ {{diag}\left( Z_{X} \right)} & 0 & 0 & 0 & 0 & {{diag}(X)} & 0 \\ 0 & {{diag}\left( Z_{S} \right)} & 0 & 0 & 0 & 0 & {{diag}(S)} \end{bmatrix} \cdot \begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \\ {\Delta\; Z_{x}} \\ {\Delta\; Z_{s}} \end{bmatrix}} = {- \begin{bmatrix} {{\nabla{f(X)}} - Z_{x} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \\ {{Z_{x} \cdot X} - {\mu\; e}} \\ {{Z_{s} \cdot S} - {\mu\; e}} \end{bmatrix}}},$ wherein [ΔX, ΔS, ΔY, Δy₁, Δy₂, ΔZ_(x), ΔZ_(s)] are unknowns, I represents an identity matrix, $\mspace{20mu}{{H_{11} = {{diag}\left( {V^{x},W^{x},V^{u},W^{U}} \right)}},\mspace{20mu}{H_{41} = \left\lbrack {{- I},R,0,0} \right\rbrack},{R = \begin{pmatrix} \left\lbrack {R_{11}I\mspace{14mu}\ldots\mspace{14mu} R_{1N_{1}}I} \right\rbrack & 0 & \ldots & 0 \\ 0 & \left\lbrack {R_{21}I\mspace{14mu}\ldots\mspace{14mu} R_{2N_{2}}I} \right\rbrack & \ldots & 0 \\ 0 & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & \left\lbrack {R_{N_{1}1}I\mspace{14mu}\ldots\mspace{14mu} R_{N_{1}N_{1}}I} \right\rbrack \end{pmatrix}},}$ a row of branch (i, j) of H₅₁ is ${\left( H_{51} \right)_{ij} = \left\lbrack {0,{- I},{\ldots\mspace{14mu} R_{ij}A_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}B_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}C_{ij}^{0}},{\ldots\mspace{14mu} R_{ij}D_{ij}^{0}},{\ldots\mspace{14mu} 0}} \right\rbrack},{X = \left( {\left\{ x_{il} \right\},\left\{ x_{ijl} \right\},\left\{ u_{il} \right\},\left\{ u_{ijl} \right\}} \right)^{T}},{{{slack}\mspace{14mu}{variables}\mspace{14mu} S} = \left( {\left\{ s_{xi} \right\},\left\{ s_{xij} \right\},\left\{ s_{ui} \right\},\left\{ s_{uij} \right\}} \right)^{T}},{{o\left( {x_{ij},x_{i}} \right)} = {x_{i\; 1} - {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{{ij}\; 1}}} - b_{x}}},{{v(X)} = {x_{{ij}\; 1} - {R_{ij} \cdot \left( {{A_{ij}^{0} \cdot u_{i\; 1}} + {B_{ij}^{0} \cdot u_{j\; 1}} + {C_{ij}^{0} \cdot u_{{ij}\; 1}} + {D_{{ij}\;}^{0} \cdot u_{{ji}\; 1}}} \right)} - b_{v}}},$ Y, y₁, and y₂ are dual variables, and Z_(x)=μ/X, Z_(s)=μ/S wherein μ is a barrier parameter.
 15. The computer readable program storage device of claim 14, wherein reducing a dimension of the linear equations further comprises solving for ΔZ _(X)=diag(X)⁻¹(−r _(Z) _(X) −diag(Z _(X))ΔX), ΔZ _(S)=diag(S)⁻¹(−r _(Z) _(S) −diag(Z _(S))ΔS), wherein r_(Z) _(X) =Z_(X)·X−μe, r_(Z) _(S) =Z_(S)·S−μe.
 16. The computer readable program storage device of claim 15, wherein the reduced dimension linear equations re-organized into primal and dual parts are $\begin{bmatrix} H_{41} & 0 & | & 0 & 0 & 0 \\ H_{51} & 0 & | & 0 & 0 & 0 \\ {- I} & {- I} & | & 0 & 0 & 0 \\  - & - & \; & - & - & - \\ {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{x} \right)}}} & 0 & | & {- I} & H_{41}^{T} & H_{51}^{T} \\ 0 & {{{diag}(S)}^{- 1}{{diag}\left( Z_{s} \right)}} & | & {- I} & 0 & 0 \end{bmatrix}{\quad{{\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \\  - \\ {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix} = \begin{bmatrix} {- r_{{y\;}_{1}}} \\ {- r_{y_{2}}} \\ {- r_{Y}} \\  - \\ {{- r_{X}} - {{{diag}(X)}^{- 1}r_{Z_{x}}}} \\ {{- r_{S}} - {{{diag}(S)}^{- 1}r_{Z_{s}}}} \end{bmatrix}},\mspace{20mu}{{{wherein}\mspace{20mu}\begin{bmatrix} r_{X} \\ r_{S} \\ r_{Y} \\ r_{y_{1}} \\ r_{y_{2}} \end{bmatrix}} = {\begin{bmatrix} {{\nabla{f(X)}} - Z_{x} - Y - {y_{1}^{T}{\nabla{o\left( {x_{ij},x_{i}} \right)}}} - {y_{2}^{T}{\nabla{v(X)}}}} \\ {{- Y} - Z_{s}} \\ {- \left( {X + S - X^{UB}} \right)} \\ {- {o\left( {x_{ij},x_{i}} \right)}} \\ {- {v(X)}} \end{bmatrix}.}}}}$
 17. The computer readable program storage device of claim 16, wherein the two systems of block matrix equations are ${{\begin{bmatrix} H_{41} & 0 \\ H_{51} & 0 \\ {- I} & {- I} \end{bmatrix}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} = \begin{bmatrix} {- r_{y\; 1}} \\ {- r_{y\; 2}} \\ {- r_{Y}} \end{bmatrix}},{{{{{and}\mspace{14mu}\begin{bmatrix} {\overset{\sim}{H}}_{11} & 0 \\ 0 & {\overset{\sim}{H}}_{22} \end{bmatrix}}\begin{bmatrix} {\Delta\; X} \\ {\Delta\; S} \end{bmatrix}} + {\begin{bmatrix} {- I} & H_{41}^{T} & H_{51}^{T} \\ {- I} & 0 & 0 \end{bmatrix}\begin{bmatrix} {\Delta\; Y} \\ {\Delta\; y_{1}} \\ {\Delta\; y_{2}} \end{bmatrix}}} = \begin{bmatrix} {- {\overset{\sim}{r}}_{X}} \\ {- {\overset{\sim}{r}}_{S}} \end{bmatrix}},{wherein}$ ${{\overset{\sim}{H}}_{11} = {H_{11} + {{{diag}(X)}^{- 1}{{diag}\left( Z_{X} \right)}}}},{{\overset{\sim}{H}}_{22} = {{{diag}(S)}^{- 1}{{diag}\left( Z_{S} \right)}}},{{\overset{\sim}{r}}_{X} = {r_{X} + {{{diag}(X)}^{- 1}r_{Z_{X}}}}},{{\overset{\sim}{r}}_{S} = {r_{S} + {{{diag}(S)}^{- 1}r_{Z_{S}}}}},{{\Delta\; X} = \begin{bmatrix} \left\{ {\Delta\; x_{i}} \right\} \\ \left\{ {\Delta\; x_{ij}} \right\} \\ \left\{ {\Delta\; u_{i}} \right\} \\ \left\{ {\Delta\; u_{ij}} \right\} \end{bmatrix}},{{\Delta\; y_{1}} = \left\{ {\Delta\; y_{1}^{x_{i}}} \right\}},{{\Delta\; y_{2}} = {\left\{ {\Delta\; y_{2}^{x_{ij}}} \right\}.}}$
 18. The computer readable program storage device of claim 17, wherein the first block matrix equation is expanded for each bus i and branch(i, j) as follows: ${{\Delta\; x_{i}} - {\sum\limits_{j}{R_{ij}\Delta\; x_{ij}}}} = r_{y_{i}}^{x_{i}}$ Δ x_(ij) − R_(ij)(A_(ij)Δ u_(i) + B_(ij)Δ u_(j) + C_(ij)Δ u_(ij) + D_(ij)Δ u_(ji)) = r_(y₂ )^(x_(ij)) Δ X + Δ S = r_(Y), and the second block matrix equation is expanded for each bus i and branch(i, j) as follows: ${{{\overset{\sim}{V}}_{i}^{X}\Delta\; x_{i}} - {\Delta\; Y^{x_{i}}} - {\Delta\; y_{1}^{x_{i}}}} = {- {\overset{\sim}{r}}_{x_{i}}}$ ${{{\overset{\sim}{W}}_{ij}^{X}\Delta\; x_{ij}} - {\Delta\; Y^{x_{ij}}} + {R_{ij}\Delta\; y_{1}^{x_{i}}} - {\Delta\; y_{2}^{x_{ij}}}} = {- {\overset{\sim}{r}}_{x_{ij}}}$ ${{{\overset{\sim}{V}}_{i}^{U}\Delta\; u_{i}} - {\Delta\; Y^{u_{i}}} + {\sum\limits_{j}{R_{ij}A_{ij}^{T}\Delta\; y_{2}^{x_{ij}}}} + {\sum\limits_{j}{R_{ji}B_{ji}^{T}\Delta\; y_{2}^{x_{ij}}}}} = {- {\overset{\sim}{r}}_{u_{i}}}$ ${{{\overset{\sim}{W}}_{ij}^{U}\Delta\; u_{ij}} - {\Delta\; Y^{u_{ij}}} + {R_{ij}C_{ij}^{T}\Delta\; y_{2}^{x_{ij}}} + {R_{ji}D_{ji}^{T}\Delta\; y_{2}^{x_{ji}}}} = {- {\overset{\sim}{r}}_{u_{ij}}}$ ${{{{\overset{\sim}{H}}_{22}\Delta\; S} - {\Delta\; Y}} = {- {\overset{\sim}{r}}_{S}}},$ and wherein the method further comprises solving for the other variables in the block matrix equations in the following order: {Δu _(i) },{Δu _(ij) }→{Δx _(ij) }→{Δx _(i) }→ΔS→ΔY→Δy ₁ →Δy ₂.
 19. The computer readable program storage device of claim 17, the method further comprising: initializing unknowns [ΔX^(k), ΔS^(k), ΔY^(k), Δy₁ ^(k), Δy₂ ^(k)] wherein k is an iteration counter; updating the barrier parameter as μ^(k+1)←ε·μ^(k), εε(0.1, . . . , 0.9); evaluating (ΔX^(k), ΔS^(k), ΔY^(k), Δy₁ ^(k), Δy₂ ^(k)) using the two block matrix equations; calculating (ΔZ_(X) ^(k),ΔZ_(S) ^(k)); determining a step length Δβ^(k); from a line search; and updating X ^(k+1) ←X ^(k)+Δβ^(k) ΔX ^(k), S ^(k+1) ←S ^(k)+Δβ^(k) ΔS ^(k), Y ^(k+1) ←Y ^(k)+Δβ^(k) ΔY ^(k), y ₁ ^(k+1) ←y ₁ ^(k)+Δβ^(k) Δy ₁ ^(k), y ₂ ^(k+1) ←y ₂ ^(k)+Δβ^(k) Δy ₂ ^(k), Z _(X) ^(k+1) ←Z _(X) ^(k)+Δβ^(k) ΔZ _(X) ^(k), Z _(X) ^(k+1) ←Z _(S) ^(k)+Δβ^(k) ΔZ _(S) ^(k), wherein the steps of updating the barrier parameter, evaluating (ΔX^(k), ΔS^(k), ΔY^(k), Δy₁ ^(k), Δy₂ ^(k)), calculating (ΔZ_(X) ^(k),ΔZ_(S) ^(k)), determining a step length Δβ^(k), and updating X^(k+1), S^(k+1), Y^(k+1), y₁ ^(k+1), y₂ ^(k+1), Z_(X) ^(k+1) and Z_(S) ^(k+1) are repeated until either a dual gap or a maximum of (r_(X), r_(S)) is less than a pre-specified minimum.
 20. The computer readable program storage device of claim 12, wherein said cost function is, ${f\left( {x,u} \right)} = {{\sum\limits_{i = 1}^{N}\left\lbrack {{\frac{1}{2}x_{i}^{T}V_{i}^{X}x_{i}} + {\frac{1}{2}{\sum\limits_{Ki}{u_{i}^{kT}V_{i}^{U}u_{i}^{k}}}} + {S_{i}^{X}x_{i}} + {S_{i}^{U}{\sum\limits_{Ki}u_{i}^{k}}}} \right\rbrack} + {\sum\limits_{i,{j = 1}}^{N}{R_{ij} \cdot \left\lbrack {{\frac{1}{2}x_{ij}^{T}W_{ij}^{X}x_{ij}} + {\frac{1}{2}{\sum\limits_{Kij}{u_{ij}^{kT}W_{ij}^{U}u_{ij}^{k}}}} + {T_{ij}^{X}x_{ij}} + {T_{ij}^{U}{\sum\limits_{Kik}u_{ij}^{k}}}} \right\rbrack}}}$ wherein bus control vector u_(i) includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(i•), bus dependent variables vector x_(i) includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector u_(ij) includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector x_(ij) includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, and connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ wherein ${x_{i} = {\sum\limits_{j = 1}^{N}{R_{ij} \cdot x_{ij}}}},$ wherein V^(X) is a Hessian of the objective function with respect to x_(i), V^(U) is a Hessian of the objective function with respect to u_(i), S^(X) is a gradient of the objective function with respect to x_(i), S^(U) is a gradient of the objective function with respect to u_(i), W^(X) is a Hessian of the objective function with respect to x_(ij), W^(U) is a Hessian of the objective function with respect to u_(ij), T^(X) is a gradient of the objective function with respect to x_(ij), and T^(U) is a gradient of the objective function with respect to u_(ij), ${subject}\mspace{14mu}{to}\mspace{14mu}{constraints}\mspace{14mu}\begin{matrix} {{x_{i}^{LB} \leq x_{i} \leq x_{i}^{UB}},} \\ {{x_{ij}^{LB} \leq x_{ij} \leq x_{ij}^{UB}},} \end{matrix}\begin{matrix} {{u_{i}^{LB} \leq {\sum\limits_{Ki}u_{i}^{k}} \leq u_{i}^{UB}},} \\ {{u_{ij}^{LB} \leq {\sum\limits_{Kij}u_{ij}^{k}} \leq u_{ij}^{UB}},} \end{matrix}$ wherein a control grid is discretized by binary variables zε{0,1} as Z _(i) ^(k) u _(ik) ^(LB) ≦u _(i) ^(k) ≦z _(i) ^(k) u _(ik) ^(UB), z _(ij) ^(k) u _(ijk) ^(LB) ≦u _(ij) ^(k) ≦z _(ij) ^(k) u _(ijk) ^(UB), wherein the binary variables are subject to constraints $z_{i}^{LB} \leq {\sum\limits_{Ki}z_{i}^{k}} \leq z_{i}^{UB}$ $z_{i}^{LB} \leq {\sum\limits_{Ki}z_{i}^{k}} \leq z_{i}^{UB}$ and satisfying a voltage law x_(ij)=R_(ij)·(Σ_(Ki)A_(ij) ^(k)u_(i) ^(k)+Σ_(Kj)B_(ij) ^(k)u_(j) ^(k)+Σ_(Kij)C_(ij) ^(k)u_(ij) ^(k)+Σ_(Kji)D_(ij) ^(k)u_(ji) ^(k)), where the matrices A_(ij) ⁰,B_(ij) ⁰,C_(ij) ⁰,D_(ij) ⁰,E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰,u_(j) ⁰,u_(ij) ⁰).
 21. The computer readable program storage device of claim 12, wherein said cost function is ƒ({tilde over (x)},u,y,v)=½[{tilde over (x)}^(T) V ^(X) {tilde over (x)}+u ^(T) V ^(U) u+y ^(T) W ^(Y) y+v ^(T) W ^(V) v]+S ^(X) {tilde over (x)}+S ^(U) u+T ^(Y) y+T ^(V) v, wherein bus control vector u includes voltage magnitudes V_(i), phase angles θ_(i), and shunt switch positions Q_(1•), bus dependent variables vector x includes net active x_(i) ^(P) and reactive x_(i) ^(Q) power injection, branch control vector v includes transformer tap selections t_(ij) and phase shifts φ_(ij), and dependent branch variable vector y includes net active x_(ij) ^(P) and reactive x_(ij) ^(Q) power injections, wherein bus dependent vector is defined as ${x = {\left. {\overset{\sim}{A}\overset{\sim}{x}}\Leftrightarrow x_{i} \right. = {{\sum\limits_{k = 1}^{G_{i}}x_{i,k}^{g}} - {\sum\limits_{k = 1}^{L_{i}}x_{i,k}^{l}} + {\sum\limits_{k = 1}^{{SH}_{i}}x_{i,k}^{sh}}}}},$ wherein where G_(i), L_(i) and SH_(i) are the number of generators, loads and shunts placed on bus i and x_(i,k) ^(g), x_(i,k) ^(l) and x_(i,k) ^(sh) are the k^(th) active and reactive generation, load and shunt variables at bus i, respectively, and Ãε{−1,0,1}^(N×2(G+L+SH)), wherein V^(X) is a Hessian of the objective function with respect to x, V^(U) is a Hessian of the objective function with respect to u, S^(X) is a gradient of the objective function with respect to x, S^(U) is a gradient of the objective function with respect to u, W^(Y) is a Hessian of the objective function with respect to y, W^(V) is a Hessian of the objective function with respect to v, T^(Y) is a gradient of the objective function with respect to y, and T^(V) is a gradient of the objective function with respect to v, said cost function having connection matrix $R_{ij} = \left\{ \begin{matrix} {1,} & {{\exists\left( {i,j} \right)},} \\ {0,} & {{o/w},} \end{matrix} \right.$ wherein ${a_{i}\overset{\sim}{x}} = {\sum\limits_{j = 1}^{N}{R_{ij}y_{ij}}}$ wherein a_(i) is i^(th) row of matrix Ã, subject to the constraints {tilde over (x)} ^(LB) ≦{tilde over (x)}≦{tilde over (x)} ^(UB) , y ^(LB) ≦y≦y ^(UB), u ^(LB) ≦u≦u ^(UB) , v ^(LB) ≦v≦v ^(UB), and satisfying a voltage law y_(ij)=R_(ij)·(A_(ij) ⁰·u_(i)+B_(ij) ⁰·u_(j)+C_(ij) ⁰·v_(ij)+D_(ij) ⁰·v_(ji)+E_(ij) ⁰), wherein the matrices A_(ij) ⁰,B_(ij) ⁰,C_(ij) ⁰,D_(ij) ⁰,E_(ij) ⁰ can be explicitly pre-calculated for every link (i, j) for a given anchor point (u_(i) ⁰,u_(j) ⁰,u_(ij) ⁰). 